diff --git a/CHANGELOG.md b/CHANGELOG.md index d200572..97cc43e 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Project.toml b/Project.toml index 7294e7a..2635914 100644 --- a/Project.toml +++ b/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 ", "Patrick Jaap "] [deps] diff --git a/examples/Example200_LowLevelPoisson.jl b/examples/Example200_LowLevelPoisson.jl index 7d568d0..7d06074 100644 --- a/examples/Example200_LowLevelPoisson.jl +++ b/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_∇) diff --git a/examples/Example205_LowLevelSpaceTimePoisson.jl b/examples/Example205_LowLevelSpaceTimePoisson.jl index 9d350f7..a6de045 100644 --- a/examples/Example205_LowLevelSpaceTimePoisson.jl +++ b/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 diff --git a/examples/Example210_LowLevelNavierStokes.jl b/examples/Example210_LowLevelNavierStokes.jl index 160b778..c626a89 100644 --- a/examples/Example210_LowLevelNavierStokes.jl +++ b/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) diff --git a/src/fedefs/h1v_br.jl b/src/fedefs/h1v_br.jl index 9a60c2d..6c8ea63 100644 --- a/src/fedefs/h1v_br.jl +++ b/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] diff --git a/src/fedefs/hdiv_bdm1.jl b/src/fedefs/hdiv_bdm1.jl index 87ea2dc..147d7d5 100644 --- a/src/fedefs/hdiv_bdm1.jl +++ b/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 diff --git a/src/fedefs/hdiv_bdm2.jl b/src/fedefs/hdiv_bdm2.jl index 1f2c664..871c42c 100644 --- a/src/fedefs/hdiv_bdm2.jl +++ b/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 diff --git a/src/fedefs/hdiv_rtk.jl b/src/fedefs/hdiv_rtk.jl index c84a894..d9ee934 100644 --- a/src/fedefs/hdiv_rtk.jl +++ b/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 diff --git a/src/feevaluator.jl b/src/feevaluator.jl index eb2179b..7191e4b 100644 --- a/src/feevaluator.jl +++ b/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, diff --git a/src/feevaluator_h1.jl b/src/feevaluator_h1.jl index ae6c7be..963dd1d 100644 --- a/src/feevaluator_h1.jl +++ b/src/feevaluator_h1.jl @@ -4,7 +4,6 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden if FEBE.subset_handler != NothingFunction subset = _update_subset!(FEBE) cvals = FEBE.cvals - fill!(cvals, 0) refbasisvals = FEBE.refbasisvals for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), k in 1:size(cvals, 1) cvals[k, dof_i, i] = refbasisvals[i][subset[dof_i], k] @@ -19,7 +18,6 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden subset = _update_subset!(FEBE) cvals = FEBE.cvals refbasisvals = FEBE.refbasisvals - fill!(cvals, 0) for i in 1:size(cvals, 3) for dof_i in 1:size(cvals, 2) cvals[1, dof_i, i] = refbasisvals[i][subset[dof_i], c] @@ -35,7 +33,6 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden coefficients = _update_coefficients!(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:size(cvals, 1) cvals[k, dof_i, i] = refbasisvals[i][subset[dof_i], k] * coefficients[k, dof_i] end @@ -48,7 +45,6 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden coefficients = _update_coefficients!(FEBE) cvals = FEBE.cvals refbasisvals = FEBE.refbasisvals - fill!(cvals, 0) for i in 1:size(cvals, 3) for dof_i in 1:size(cvals, 2) cvals[1, dof_i, i] = refbasisvals[i][subset[dof_i], c] * coefficients[c, dof_i] @@ -65,10 +61,9 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - fill!(cvals, 0) - for i in 1:size(cvals, 3), c in 1:length(offsets), j in 1:size(L2GAinv, 1), k in 1:size(L2GAinv, 2), dof_i in 1:size(cvals, 2) + for i in 1:size(cvals, 3), c in 1:length(offsets), k in 1:size(L2GAinv, 2), dof_i in 1:size(cvals, 2) # compute duc/dxk - cvals[k + offsets[c], dof_i, i] += L2GAinv[k, j] * refbasisderivvals[subset[dof_i] + offsets2[c], j, i] + cvals[k + offsets[c], dof_i, i] = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, i]) end return nothing end @@ -82,13 +77,8 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - fill!(cvals, 0) for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), c in 1:length(offsets), k in 1:size(L2GAinv, 2) - for j in 1:size(L2GAinv, 1) - # compute duc/dxk - cvals[k + offsets[c], dof_i, i] += L2GAinv[k, j] * refbasisderivvals[subset[dof_i] + offsets2[c], j, i] - end - cvals[k + offsets[c], dof_i, i] *= coefficients[c, dof_i] + cvals[k + offsets[c], dof_i, i] = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, i]) * coefficients[c, dof_i] end return nothing end @@ -103,15 +93,15 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Symm offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals fill!(cvals, 0) - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), c in 1:length(offsets), k in 1:size(L2GAinv, 2), j in 1:size(L2GAinv, 1) + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), c in 1:length(offsets), k in 1:size(L2GAinv, 2) # compute duc/dxk and put it into the right spot in the Voigt vector + temp = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, i]) if k != c - cvals[compression[k + offsets[c]], dof_i, i] += offdiagval * L2GAinv[k, j] * refbasisderivvals[subset[dof_i] + offsets2[c], j, i] + cvals[compression[k + offsets[c]], dof_i, i] += offdiagval * temp else - cvals[compression[k + offsets[c]], dof_i, i] += L2GAinv[k, j] * refbasisderivvals[subset[dof_i] + offsets2[c], j, i] + cvals[compression[k + offsets[c]], dof_i, i] += temp end end - return nothing end @@ -123,8 +113,8 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Dive offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals fill!(cvals, 0) - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), k in 1:size(L2GAinv, 2), j in 1:size(L2GAinv, 1) - cvals[1, dof_i, i] += L2GAinv[k, j] * refbasisderivvals[subset[dof_i] + offsets2[k], j, i] + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), k in 1:size(L2GAinv, 2) + cvals[1, dof_i, i] += @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[k], :, i]) end return nothing end @@ -138,8 +128,8 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Dive offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals fill!(cvals, 0) - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), k in 1:size(L2GAinv, 2), j in 1:size(L2GAinv, 1) - cvals[1, dof_i, i] += L2GAinv[k, j] * refbasisderivvals[subset[dof_i] + offsets2[k], j, i] * coefficients[k, dof_i] + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), k in 1:size(L2GAinv, 2) + cvals[1, dof_i, i] += @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[k], :, i]) * coefficients[k, dof_i] end return nothing end @@ -167,7 +157,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Tang 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) + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) cvals[1, dof_i, i] = refbasisvals[i][subset[dof_i], 1] * normal[2] - refbasisvals[i][subset[dof_i], 2] * normal[1] end return nothing @@ -288,14 +278,9 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - fill!(cvals, 0) - for i in 1:size(cvals, 3) - for dof_i in 1:size(cvals, 2) - for j in 1:size(L2GAinv, 2) - cvals[1, dof_i, i] -= L2GAinv[2, j] * refbasisderivvals[subset[dof_i], j, i] # -du/dy - cvals[2, dof_i, i] += L2GAinv[1, j] * refbasisderivvals[subset[dof_i], j, i] # du/dx - end - end + @views for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] = -dot(L2GAinv[2, :], refbasisderivvals[subset[dof_i], :, i]) # -du1/dy + cvals[2, dof_i, i] = dot(L2GAinv[1, :], refbasisderivvals[subset[dof_i], :, i]) # du2/dx end return nothing end @@ -308,14 +293,9 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - fill!(cvals, 0) - for i in 1:size(cvals, 3) - for dof_i in 1:size(cvals, 2) - for j in 1:size(L2GAinv, 2) - cvals[1, dof_i, i] -= L2GAinv[2, j] * refbasisderivvals[subset[dof_i], j, i] # -du1/dy - cvals[1, dof_i, i] += L2GAinv[1, j] * refbasisderivvals[subset[dof_i] + offsets2[2], j, i] # du2/dx - end - end + @views for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] = dot(L2GAinv[2, :], refbasisderivvals[subset[dof_i] + offsets2[1], :, i]) # -du1/dy + cvals[1, dof_i, i] += dot(L2GAinv[1, :], refbasisderivvals[subset[dof_i] + offsets2[2], :, i]) # du2/dx end return nothing end @@ -329,18 +309,13 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - fill!(cvals, 0) - for i in 1:size(cvals, 3) - for dof_i in 1:size(cvals, 2) - for k in 1:3 - cvals[1, dof_i, i] += L2GAinv[2, k] * refbasisderivvals[subset[dof_i] + offsets2[3], k, i] # du3/dx2 - cvals[1, dof_i, i] -= L2GAinv[3, k] * refbasisderivvals[subset[dof_i] + offsets2[2], k, i] # - du2/dx3 - cvals[2, dof_i, i] += L2GAinv[3, k] * refbasisderivvals[subset[dof_i] + offsets2[1], k, i] # du1/dx3 - cvals[2, dof_i, i] -= L2GAinv[1, k] * refbasisderivvals[subset[dof_i] + offsets2[3], k, i] # - du3/dx1 - cvals[3, dof_i, i] += L2GAinv[1, k] * refbasisderivvals[subset[dof_i] + offsets2[2], k, i] # du2/dx1 - cvals[3, dof_i, i] -= L2GAinv[2, k] * refbasisderivvals[subset[dof_i] + offsets2[1], k, i] # - du1/dx2 - end - end + @views for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] = dot(L2GAinv[2, :], refbasisderivvals[subset[dof_i] + offsets2[3], :, i]) # du3/dx2 + cvals[1, dof_i, i] -= dot(L2GAinv[3, :], refbasisderivvals[subset[dof_i] + offsets2[2], :, i]) # - du2/dx3 + cvals[2, dof_i, i] = dot(L2GAinv[3, :], refbasisderivvals[subset[dof_i] + offsets2[1], :, i]) # du1/dx3 + cvals[2, dof_i, i] -= dot(L2GAinv[1, :], refbasisderivvals[subset[dof_i] + offsets2[3], :, i]) # - du3/dx1 + cvals[3, dof_i, i] = dot(L2GAinv[1, :], refbasisderivvals[subset[dof_i] + offsets2[2], :, i]) # du2/dx1 + cvals[3, dof_i, i] -= dot(L2GAinv[2, :], refbasisderivvals[subset[dof_i] + offsets2[1], :, i]) # - du1/dx2 end return nothing end @@ -363,10 +338,8 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Tang for i in 1:size(cvals, 3) for dof_i in 1:size(cvals, 2) for c in 1:length(offsets), k in 1:size(L2GAinv, 1) - for j in 1:size(L2GAinv, 2) - # compute duc/dxk - cvals[1, dof_i, i] += L2GAinv[k, j] * refbasisderivvals[subset[dof_i] + offsets2[c], j, i] * tangent[c] - end + temp = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, i]) + cvals[1, dof_i, i] += temp * tangent[c] end end end diff --git a/src/feevaluator_hcurl.jl b/src/feevaluator_hcurl.jl index f5a6806..67ea42b 100644 --- a/src/feevaluator_hcurl.jl +++ b/src/feevaluator_hcurl.jl @@ -5,12 +5,8 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden coefficients = _update_coefficients!(FEBE) refbasisvals = FEBE.refbasisvals cvals = FEBE.cvals - fill!(cvals, 0) for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2), k in 1:size(cvals, 1) - for l in 1:size(L2GAinv, 2) - cvals[k, dof_i, i] += L2GAinv[k, l] * refbasisvals[i][subset[dof_i], l] - end - cvals[k, dof_i, i] *= coefficients[k, dof_i] + cvals[k, dof_i, i] = @views dot(L2GAinv[k, :], refbasisvals[i][subset[dof_i], :]) * coefficients[k, dof_i] end return nothing end @@ -43,8 +39,8 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl fill!(cvals, 0) det = FEBE.L2G.det for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - cvals[1, dof_i, i] = refbasisderivvals[subset[dof_i] + offsets2[1], 2, i] - cvals[1, dof_i, i] -= refbasisderivvals[subset[dof_i] + offsets2[2], 1, i] + cvals[1, dof_i, i] = refbasisderivvals[subset[dof_i] + offsets2[2], 1, i] # du2/dx + cvals[1, dof_i, i] -= refbasisderivvals[subset[dof_i] + offsets2[1], 2, i] # -du1/dy cvals[1, dof_i, i] *= coefficients[1, dof_i] / det end return nothing diff --git a/src/feevaluator_hdiv.jl b/src/feevaluator_hdiv.jl index 6419842..62ac927 100644 --- a/src/feevaluator_hdiv.jl +++ b/src/feevaluator_hdiv.jl @@ -6,13 +6,9 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden det = FEBE.L2G.det # 1 alloc cvals = FEBE.cvals refbasisvals = FEBE.refbasisvals - fill!(cvals, 0) for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) for k in 1:size(L2GM, 1) - for l in 1:size(L2GM, 2) - cvals[k, dof_i, i] += L2GM[k, l] * refbasisvals[i][subset[dof_i], l] - end - cvals[k, dof_i, i] *= coefficients[k, dof_i] / det + cvals[k, dof_i, i] = @views dot(L2GM[k, :], refbasisvals[i][subset[dof_i], :]) * coefficients[k, dof_i] / det end end return nothing @@ -27,12 +23,8 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden det = FEBE.L2G.det # 1 alloc cvals = FEBE.cvals refbasisvals = FEBE.refbasisvals - fill!(cvals, 0) for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - for l in 1:size(L2GM, 2) - cvals[1, dof_i, i] += L2GM[c, l] * refbasisvals[i][subset[dof_i], l] - end - cvals[1, dof_i, i] *= coefficients[c, dof_i] / det + cvals[1, dof_i, i] = @views dot(L2GM[c, :], refbasisvals[i][subset[dof_i], :]) * coefficients[c, dof_i] / det end return nothing end @@ -96,12 +88,19 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad fill!(cvals, 0) det = FEBE.L2G.det # 1 alloc for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - for c in 1:size(L2GM, 1), k in 1:size(L2GAinv, 1) - # compute duc/dxk - for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2) - cvals[k + offsets[c], dof_i, i] += L2GAinv[k, m] * L2GM[c, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] + for k in 1:size(L2GAinv, 1) + for j in 1:size(L2GM, 2) + # apply inverse transform to d/dx_k of j-th component of i-th reference basis function + temp = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[j], :, i]) + # add contribution d/dx_k of c-th component of i-th reference basis function + for c in 1:size(L2GM, 1) + cvals[k + offsets[c], dof_i, i] += L2GM[c, j] * temp + end + end + # apply trafo and orientation factors + for c in 1:size(L2GM, 1) + cvals[k + offsets[c], dof_i, i] *= coefficients[c, dof_i] / det end - cvals[k + offsets[c], dof_i, i] *= coefficients[c, dof_i] / det end end return nothing @@ -119,10 +118,12 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl refbasisderivvals = FEBE.refbasisderivvals fill!(cvals, 0) det = FEBE.L2G.det # 1 alloc - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2) - cvals[1, dof_i, i] -= L2GAinv[2, m] * L2GM[1, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] / det - cvals[1, dof_i, i] += L2GAinv[1, m] * L2GM[2, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det + for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2) + A = L2GAinv[2, m] * L2GM[1, j] / det + B = L2GAinv[1, m] * L2GM[2, j] / det + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] -= A * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] + cvals[1, dof_i, i] += B * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] end end return nothing @@ -140,14 +141,20 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl refbasisderivvals = FEBE.refbasisderivvals fill!(cvals, 0) det = FEBE.L2G.det # 1 alloc - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2) - cvals[1, dof_i, i] -= L2GAinv[3, m] * L2GM[2, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det - cvals[1, dof_i, i] += L2GAinv[2, m] * L2GM[3, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[3, dof_i] / det - cvals[2, dof_i, i] -= L2GAinv[1, m] * L2GM[3, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[3, dof_i] / det - cvals[2, dof_i, i] += L2GAinv[3, m] * L2GM[1, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] / det - cvals[3, dof_i, i] -= L2GAinv[2, m] * L2GM[1, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] / det - cvals[3, dof_i, i] += L2GAinv[1, m] * L2GM[2, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det + for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2) + A = L2GAinv[3, m] * L2GM[2, j] / det + B = L2GAinv[2, m] * L2GM[3, j] / det + C = L2GAinv[1, m] * L2GM[3, j] / det + D = L2GAinv[3, m] * L2GM[1, j] / det + E = L2GAinv[2, m] * L2GM[1, j] / det + F = L2GAinv[1, m] * L2GM[2, j] / det + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] -= A * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] + cvals[1, dof_i, i] += B * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[3, dof_i] + cvals[2, dof_i, i] -= C * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[3, dof_i] + cvals[2, dof_i, i] += D * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] + cvals[3, dof_i, i] -= E * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] + cvals[3, dof_i, i] += F * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] end end return nothing diff --git a/src/interpolators.jl b/src/interpolators.jl index 4c5174f..b7ba21b 100644 --- a/src/interpolators.jl +++ b/src/interpolators.jl @@ -392,7 +392,7 @@ function MomentInterpolator( end weights, xref = QF.w, QF.xref nweights = length(weights) - for item::Int in items + for item::Ti in items if item < 1 continue end @@ -609,7 +609,7 @@ function FunctionalInterpolator( end weights, xref = QF.w, QF.xref nweights = length(weights) - for item::Int in items + for item::Ti in items if item < 1 continue end diff --git a/src/lazy_interpolate.jl b/src/lazy_interpolate.jl index 24b22fc..65ba1fb 100644 --- a/src/lazy_interpolate.jl +++ b/src/lazy_interpolate.jl @@ -62,10 +62,10 @@ function lazy_interpolate!( operators = [(1, Identity)]; postprocess = standard_kernel, xtrafo = nothing, - items = [], + items = Ti[], resultdim = get_ncomponents(eltype(target.FES)), not_in_domain_value = 1.0e30, - start_cell = 1, + start_cell::Ti = Ti(1), only_localsearch = false, use_cellparents::Bool = false, eps = 1.0e-13, @@ -82,8 +82,8 @@ function lazy_interpolate!( PE = PointEvaluator(postprocess, operators, source; TCoeff = T1) xref = zeros(Tv, xdim_source) x_source = zeros(Tv, xdim_source) - cell::Int = start_cell - lastnonzerocell::Int = start_cell + cell::Ti = start_cell + lastnonzerocell::Ti = start_cell same_cells::Bool = xgrid == target.FES.xgrid CF::CellFinder{Tv, Ti} = CellFinder(xgrid) diff --git a/src/plots.jl b/src/plots.jl index f9f5f41..30a4e91 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -5,8 +5,8 @@ A standard scalarplot of (the operator evaluation of) a finite element vector in All kwargs of the calling method are transferred to the scalarplot. """ -function GridVisualize.scalarplot!(vis::Union{Nothing, Dict{Symbol, Any}}, feVectorBlock::FEVectorBlock, operator = Identity; abs = false, component = 1, kwargs...) - return GridVisualize.scalarplot!(vis, feVectorBlock.FES.dofgrid, view(nodevalues(feVectorBlock, operator; abs = abs), component, :); kwargs...) +function GridVisualize.scalarplot!(vis::Union{Nothing, Dict{Symbol, Any}}, feVectorBlock::FEVectorBlock, operator = Identity; abs = false, component = 1, title = feVectorBlock.name, kwargs...) + return GridVisualize.scalarplot!(vis, feVectorBlock.FES.dofgrid, view(nodevalues(feVectorBlock, operator; abs = abs), component, :); title = title, kwargs...) end """ scalarplot(feVectorBlock::FEVectorBlock, operator = Identity; kwargs...) @@ -29,7 +29,7 @@ independently. Thus, a discontinuous plot is generated. All kwargs of the calling method are transferred to the scalarplot. """ -function broken_scalarplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...) +function broken_scalarplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name, kwargs...) dofgrid = feVectorBlock.FES.dofgrid cell_nodes = dofgrid[CellNodes] @@ -73,8 +73,8 @@ A standard vectorplot of (the operator evaluation of) a finite element vector. All kwargs of the calling method are transferred to the vectorplot. """ -function GridVisualize.vectorplot(feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name, kwargs...) - return GridVisualize.vectorplot(feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); title = title, kwargs...) +function GridVisualize.vectorplot(feVectorBlock::FEVectorBlock, operator = Identity; kwargs...) + return GridVisualize.vectorplot(feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); kwargs...) end @@ -88,6 +88,6 @@ All kwargs of the calling method are transferred to the streamplot. function GridVisualize.streamplot!(p, feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name * " (streamlines)", kwargs...) return streamplot!(p, feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); title = title, kwargs...) end -function GridVisualize.streamplot(feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name * " (streamlines)", kwargs...) - return streamplot(feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); title = title, kwargs...) +function GridVisualize.streamplot(feVectorBlock::FEVectorBlock, operator = Identity; kwargs...) + return streamplot(feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); kwargs...) end diff --git a/src/point_evaluator.jl b/src/point_evaluator.jl index 2aa5743..852ecfa 100644 --- a/src/point_evaluator.jl +++ b/src/point_evaluator.jl @@ -1,11 +1,11 @@ -mutable struct PointEvaluator{Tv <: Real, TCoeff <: Real, UT, KFT <: Function} +mutable struct PointEvaluator{Tv <: Real, Ti <: Integer, TCoeff <: Real, UT, KFT <: Function} u_args::Array{UT, 1} ops_args::Array{DataType, 1} kernel::KFT BE_args::Any L2G::Any CF::Any - lastitem::Int + lastitem::Ti eval_selector::Any evaluator_bary::Any evaluator::Any @@ -61,11 +61,11 @@ $(_myprint(default_peval_kwargs())) After construction, call `initialize!` to prepare the evaluator for a given solution, then use `evaluate!` or `evaluate_bary!` to perform point evaluations. """ -function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, TCoeff = Float64, kwargs...) +function PointEvaluator(kernel, u_args, ops_args, sol = nothing; Tv = Float64, Ti = Int32, TCoeff = Float64, kwargs...) parameters = Dict{Symbol, Any}(k => v[1] for (k, v) in default_peval_kwargs()) _update_params!(parameters, kwargs) @assert length(u_args) == length(ops_args) - PE = PointEvaluator{Tv, TCoeff, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, 1, nothing, nothing, nothing, zeros(Tv, 2), parameters) + PE = PointEvaluator{Tv, Ti, TCoeff, typeof(u_args[1]), typeof(kernel)}(u_args, ops_args, kernel, nothing, nothing, nothing, Ti(1), nothing, nothing, nothing, zeros(Tv, 2), parameters) if sol !== nothing initialize!(PE, sol) end @@ -121,7 +121,7 @@ $(_myprint(default_peval_kwargs())) # Notes - This function must be called before using `evaluate!` or `evaluate_bary!` with the `PointEvaluator`. """ -function initialize!(O::PointEvaluator{T, TCoeff, UT}, sol; time = 0, kwargs...) where {T, TCoeff, UT} +function initialize!(O::PointEvaluator{T, Ti, TCoeff, UT}, sol; time = 0, kwargs...) where {T, Ti, TCoeff, UT} _update_params!(O.parameters, kwargs) if UT <: Integer ind_args = O.u_args @@ -131,7 +131,6 @@ function initialize!(O::PointEvaluator{T, TCoeff, UT}, sol; time = 0, kwargs...) FES_args = [sol[j].FES for j in ind_args] nargs = length(FES_args) xgrid = FES_args[1].xgrid - Ti = eltype(xgrid[CellNodes]) EGs = xgrid[UniqueCellGeometries] AT = ON_CELLS gridAT = EffAT4AssemblyType(get_AT(FES_args[1]), AT) @@ -172,10 +171,10 @@ function initialize!(O::PointEvaluator{T, TCoeff, UT}, sol; time = 0, kwargs...) function _evaluate_bary!( result, BE_args::Array{<:FEEvaluator, 1}, - L2G::L2GTransformer, + L2G::L2GTransformer{Tv, Ti}, xref, item, # cell used to evaluate local coordinates - ) + ) where {Tv, Ti} for id in 1:nargs # update basis evaluations at xref @@ -246,10 +245,10 @@ Evaluates the PointEvaluator at the specified reference coordinates in the cell """ function evaluate_bary!( result, - PE::PointEvaluator, + PE::PointEvaluator{Tv, Ti}, xref, item, - ) + ) where {Tv, Ti} ## find cell geometry id j = PE.eval_selector(item) diff --git a/src/qpinfos.jl b/src/qpinfos.jl index 5f57711..33f3418 100644 --- a/src/qpinfos.jl +++ b/src/qpinfos.jl @@ -9,7 +9,7 @@ A mutable struct that encapsulates information about the current quadrature poin - `region::Ti`: Index of the current region (with respect to the current assembly type). - `volume::TvG`: Volume associated with the current item. - `normal::Vector{TvG}`: Outward normal vector at the quadrature point, if applicable (e.g., for face assembly/integration). -- `time::Ttime`: Current time value, useful for time-dependent problems. +- `time::Ttime`: Current time value, useful for time-dependent problems (default: time = 1.0). - `x::Vector{Tx}`: World (physical) coordinates of the quadrature point. - `xref::Vector{Txref}`: Reference (local) coordinates of the quadrature point within the reference element. - `grid::ExtendableGrid{TvG, TiG}`: Reference to the underlying grid or mesh structure. diff --git a/src/quadrature.jl b/src/quadrature.jl index f64c720..0f4e207 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -767,7 +767,7 @@ function integrate!( resultdim::Int = (typeof(integral4items) <: AbstractArray{T, 1}) ? length(offset) : size(integral4items, 1) result::Vector{T} = zeros(T, resultdim) return if typeof(integral4items) <: AbstractArray{T, 1} - function _integrate_cell_1d!(integral4items, ilocal2global::L2GTransformer{T}, iqf::QuadratureRule{T}, QP, xCellParents, item) + function _integrate_cell_1d!(integral4items, ilocal2global::L2GTransformer{T}, iqf::QuadratureRule{T}, QP, xCellParents, item::Ti) update_trafo!(ilocal2global, item) QP.item = item @@ -787,7 +787,7 @@ function integrate!( end fill!(integral4items, 0) - for item::Int in items + for item::Ti in items if xItemRegions[item] > 0 if !(visit_region[xItemRegions[item]]) || AT == ON_IFACES continue @@ -807,7 +807,7 @@ function integrate!( _integrate_cell_1d!(integral4items, local2global[iEG], qf[iEG], QP, xCellParents, item) end else # <: AbstractArray{T,2} - function _integrate_cell_2d!(integral4items, ilocal2global::L2GTransformer{T}, iqf::QuadratureRule{T}, QP, xCellParents, item) + function _integrate_cell_2d!(integral4items, ilocal2global::L2GTransformer{T}, iqf::QuadratureRule{T}, QP, xCellParents, item::Ti) update_trafo!(ilocal2global, item) QP.item = item @@ -827,7 +827,7 @@ function integrate!( end fill!(integral4items, 0) - for item::Int in items + for item::Ti in items if xItemRegions[item] > 0 if !(visit_region[xItemRegions[item]]) || AT == ON_IFACES continue diff --git a/src/reconstructionoperators.jl b/src/reconstructionoperators.jl index 91d326c..716ee51 100644 --- a/src/reconstructionoperators.jl +++ b/src/reconstructionoperators.jl @@ -7,7 +7,7 @@ DefaultName4Operator(::Type{Reconstruct{FETypeR, O}}) where {FETypeR, O} = "R(" DefaultName4Operator(::Type{WeightedReconstruct{FETypeR, O, w}}) where {FETypeR, O, w} = "wR(r" * DefaultName4Operator(O) * ")" struct FEReconstEvaluator{T, TvG, TiG, FEType, FEType2, stdop, RH} <: 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) FEB::SingleFEEvaluator{T, TvG, TiG, stdop, FEType2} # FEBasisEvaluator for stdop in reconstruction space cvals::Array{T, 3} # current operator vals on item (reconstruction) diff --git a/test/test_operators.jl b/test/test_operators.jl index a943e8e..48c2c71 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -4,13 +4,17 @@ function run_operator_tests() println("============================") println("Testing Operator Evaluations") println("============================") - error = test_derivatives2D() + error = test_derivatives2D(H1P2{2, 2}, 2) @test error < 1.0e-14 - error = test_derivatives2D_hdiv() + error = test_derivatives2D(HDIVBDM2{2}, 2) @test error < 1.0e-14 - error = test_derivatives3D() + error = test_derivatives2D(HCURLN1{2}, 1) @test error < 1.0e-14 - error = test_derivatives3D_hdiv() + error = test_derivatives3D(H1P2{3, 3}, 2) + @test error < 1.0e-14 + error = test_derivatives3D(HDIVRT1{3}, 1) + @test error < 1.0e-14 + error = test_derivatives3D(HCURLN0{3}, 1) @test error < 1.0e-14 test_reconstructions() end @@ -64,161 +68,200 @@ function test_reconstructions() return end -function test_derivatives2D() +function test_derivatives2D(fetype, order) ## define test function and expected operator evals function testf(result, qpinfo) x = qpinfo.x - result[1] = x[1]^2 - return result[2] = 3 * x[2]^2 + x[1] * x[2] + if order == 2 + result[1] = x[1]^2 + result[2] = 3 * x[2]^2 + x[1] * x[2] + elseif order == 1 + result[1] = x[1] + x[2] + 1 + result[2] = 3 * x[2] - x[1] + end + return nothing end - ## expected values of operators in cell midpoint - expected_L = [2, 6] # expected Laplacian - expected_H = [2, 0, 0, 0, 0, 1, 1, 6] # expected Hessian - expected_symH = [2, 0, 0, 0, 6, 1] # expected symmetric Hessian - expected_symH2 = [2, 0, 0, 0, 6, sqrt(2)] # expected symmetric Hessian - expected_curl2 = [1 / 3] + ## define grid = a single non-reference triangle + xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]') # midpoint = [0.0, 1 / 3] - ## define grid = a single non-refenrece triangle - xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]') + ## expected values of operators in cell midpoint + if order == 2 + expected_id = [0, 1 / 3] + expected_L = [2, 6] # expected Laplacian + expected_H = [2, 0, 0, 0, 0, 1, 1, 6] # expected Hessian + expected_symH = [2, 0, 0, 0, 6, 1] # expected symmetric Hessian + expected_symH2 = [2, 0, 0, 0, 6, sqrt(2)] # expected symmetric Hessian + expected_curl2 = [1 / 3] + expected_grad = [0, 0, 1 / 3, 2] + expected_div = [2] + elseif order == 1 + expected_id = [4 / 3, 1] + expected_curl2 = [-2] + end ## define P2-Courant finite element space - FES = FESpace{H1P2{2, 2}}(xgrid) + FES = FESpace{fetype}(xgrid) show(devnull, FES) + ## interpolate quadratic testfunction + Iu = FEVector(FES) + interpolate!(Iu[1], testf) + ## get midpoint quadrature rule for constants qf = QuadratureRule{Float64, Triangle2D}(0) - ## define FE basis Evaluator for Hessian + FEBE_id = FEEvaluator(FES, Identity, qf) FEBE_curl2 = FEEvaluator(FES, Curl2D, qf) + update_basis!(FEBE_id, 1) + update_basis!(FEBE_curl2, 1) + ## check if operator evals have the correct length + @assert size(FEBE_id.cvals, 1) == length(expected_id) + @assert size(FEBE_curl2.cvals, 1) == length(expected_curl2) + # evaluate at quadrature points = cell midpoint + id = zeros(Float64, 2) + curl2 = zeros(Float64, 1) + eval_febe!(id, FEBE_id, Iu.entries[FES[CellDofs][:, 1]], 1) + eval_febe!(curl2, FEBE_curl2, Iu.entries[FES[CellDofs][:, 1]], 1) + ## compute errors to expected values + error_id = sqrt(sum((id - expected_id) .^ 2)) + error_curl2 = sqrt(sum((curl2 - expected_curl2) .^ 2)) + println("EG = Triangle2D | $fetype | operator = Identity | error = $error_id") + println("EG = Triangle2D | $fetype | operator = Curl2 | error = $error_curl2") + if fetype <: AbstractHcurlFiniteElement + return maximum([error_id, error_curl2]) + end + + grad = zeros(Float64, 4) + div = zeros(Float64, 1) + FEBE_grad = FEEvaluator(FES, Gradient, qf) + FEBE_div = FEEvaluator(FES, Divergence, qf) + update_basis!(FEBE_grad, 1) + update_basis!(FEBE_div, 1) + @assert size(FEBE_grad.cvals, 1) == length(expected_grad) + @assert size(FEBE_div.cvals, 1) == length(expected_div) + eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) + eval_febe!(div, FEBE_div, Iu.entries[FES[CellDofs][:, 1]], 1) + error_grad = sqrt(sum((grad - expected_grad) .^ 2)) + error_div = sqrt(sum((div - expected_div) .^ 2)) + println("EG = Triangle2D | $fetype | operator = Gradient | error = $error_grad") + println("EG = Triangle2D | $fetype | operator = Divergence | error = $error_div") + + if fetype <: AbstractHdivFiniteElement + return maximum([error_id, error_curl2, error_grad]) + end + + # do the same for second order derivatives FEBE_L = FEEvaluator(FES, Laplacian, qf) FEBE_H = FEEvaluator(FES, Hessian, qf) FEBE_symH = FEEvaluator(FES, SymmetricHessian{1}, qf) FEBE_symH2 = FEEvaluator(FES, SymmetricHessian{sqrt(2)}, qf) - - ## update on cell 1 update_basis!(FEBE_L, 1) update_basis!(FEBE_H, 1) update_basis!(FEBE_symH, 1) update_basis!(FEBE_symH2, 1) - update_basis!(FEBE_curl2, 1) - - ## interpolate quadratic testfunction - Iu = FEVector(FES) - interpolate!(Iu[1], testf) - - ## check if operator evals have the correct length @assert size(FEBE_L.cvals, 1) == length(expected_L) @assert size(FEBE_H.cvals, 1) == length(expected_H) @assert size(FEBE_symH.cvals, 1) == length(expected_symH) @assert size(FEBE_symH2.cvals, 1) == length(expected_symH2) - @assert size(FEBE_curl2.cvals, 1) == length(expected_curl2) - - ## eval 2nd order derivatives at only quadrature point 1 - ## since function is quadratic this should be constant H = zeros(Float64, 8) symH = zeros(Float64, 6) symH2 = zeros(Float64, 6) L = zeros(Float64, 2) - curl2 = zeros(Float64, 1) eval_febe!(L, FEBE_L, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(H, FEBE_H, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(symH, FEBE_symH, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(symH2, FEBE_symH2, Iu.entries[FES[CellDofs][:, 1]], 1) - eval_febe!(curl2, FEBE_curl2, Iu.entries[FES[CellDofs][:, 1]], 1) - - ## compute errors to expected values error_L = sqrt(sum((L - expected_L) .^ 2)) error_H = sqrt(sum((H - expected_H) .^ 2)) error_symH = sqrt(sum((symH - expected_symH) .^ 2)) error_symH2 = sqrt(sum((symH2 - expected_symH2) .^ 2)) - error_curl2 = sqrt(sum((curl2 - expected_curl2) .^ 2)) - println("EG = Triangle2D | operator = Curl2 | error = $error_curl2") - println("EG = Triangle2D | operator = Laplacian | error = $error_L") - println("EG = Triangle2D | operator = Hessian | error = $error_H") - println("EG = Triangle2D | operator = SymmetricHessian{1} | error = $error_symH") - println("EG = Triangle2D | operator = SymmetricHessian{√2} | error = $error_symH2") - - return maximum([error_curl2, error_L, error_H, error_symH, error_symH2]) + println("EG = Triangle2D | $fetype | operator = Laplacian | error = $error_L") + println("EG = Triangle2D | $fetype | operator = Hessian | error = $error_H") + println("EG = Triangle2D | $fetype | operator = SymmetricHessian{1} | error = $error_symH") + println("EG = Triangle2D | $fetype | operator = SymmetricHessian{√2} | error = $error_symH2") + return maximum([error_id, error_curl2, error_L, error_H, error_symH, error_symH2, error_grad, error_div]) end - -function test_derivatives2D_hdiv() +function test_derivatives3D(fetype, order) ## define test function and expected operator evals function testf(result, qpinfo) x = qpinfo.x - result[1] = x[1]^2 - return result[2] = 3 * x[2]^2 + x[1] * x[2] + if order == 2 + result[1] = x[1]^2 + x[3] * x[2] + result[2] = 3 * x[3]^2 + x[1] * x[2] + result[3] = x[1] * x[2] + elseif order == 1 + result[1] = x[2] + result[2] = 3 * x[3] + result[3] = x[2] + end + return nothing end - ## expected values of operators in cell midpoint - expected_curl2 = [1 / 3] - ## define grid = a single non-refenrece triangle - xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]') + xgrid = reference_domain(Tetrahedron3D) + xgrid[Coordinates][:, 2] = [2, 0, 0] # midpoint = [0.5, 0.25, 0.25] ## define P2-Courant finite element space - FES = FESpace{HDIVBDM2{2}}(xgrid) + FES = FESpace{fetype}(xgrid) show(devnull, FES) - ## get midpoint quadrature rule for constants - qf = QuadratureRule{Float64, Triangle2D}(0) - - ## define FE basis Evaluator for Hessian - FEBE_curl2 = FEEvaluator(FES, Curl2D, qf) - - ## update on cell 1 - update_basis!(FEBE_curl2, 1) - ## interpolate quadratic testfunction Iu = FEVector(FES) interpolate!(Iu[1], testf) - ## check if operator evals have the correct length - @assert size(FEBE_curl2.cvals, 1) == length(expected_curl2) - - # compute derivatives - curl2 = zeros(Float64, 1) - eval_febe!(curl2, FEBE_curl2, Iu.entries[FES[CellDofs][:, 1]], 1) + ## expected values of operators in cell midpoint + if order == 2 + expected_L = [2, 6, 0] # expected Laplacian + expected_H = [2, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 6, 0, 1, 0, 1, 0, 0, 0, 0, 0] # expected Hessian + expected_symH = [2, 0, 0, 1, 0, 0, 0, 0, 6, 0, 0, 1, 0, 0, 0, 0, 0, 1] # expected symmetric Hessian + expected_symH2 = [2, 0, 0, sqrt(2), 0, 0, 0, 0, 6, 0, 0, sqrt(2), 0, 0, 0, 0, 0, sqrt(2)] # expected symmetric Hessian + expected_curl3 = [1 / 2 - 6 / 4, 0, 1 / 4 - 1 / 4] + expected_grad = [1, 0.25, 0.25, 0.25, 0.5, 1.5, 0.25, 0.5, 0] # expected Gradient + expected_div = [1.5] + elseif order == 1 + expected_curl3 = [1 - 3, 0, -1] + expected_grad = [0, 1, 0, 0, 0, 3, 0, 1, 0] # expected Gradient + expected_div = [0] + end - ## compute errors to expected values - error_curl2 = sqrt(sum((curl2 - expected_curl2) .^ 2)) - println("EG = Triangle2D | operator = Curl2 | error = $error_curl2") + ## get midpoint quadrature rule for constants + qf = QuadratureRule{Float64, Tetrahedron3D}(0) - return maximum([error_curl2]) -end + FEBE_curl3 = FEEvaluator(FES, Curl3D, qf) + update_basis!(FEBE_curl3, 1) + @assert size(FEBE_curl3.cvals, 1) == length(expected_curl3) + FEBE_grad = FEEvaluator(FES, Gradient, qf) + FEBE_div = FEEvaluator(FES, Divergence, qf) + curl3 = zeros(Float64, 3) + eval_febe!(curl3, FEBE_curl3, Iu.entries[FES[CellDofs][:, 1]], 1) + error_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2)) + println("EG = Tetrahedron3D | $fetype | operator = Curl3 | error = $error_curl3") -function test_derivatives3D() - ## define test function and expected operator evals - function testf(result, qpinfo) - x = qpinfo.x - result[1] = x[1]^2 + x[3] * x[2] - result[2] = 3 * x[3]^2 + x[1] * x[2] - return result[3] = x[1] * x[2] + if fetype <: AbstractHcurlFiniteElement + return maximum([error_curl3]) end - ## expected values of operators in cell midpoint - expected_L = [2, 6, 0] # expected Laplacian - expected_H = [2, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 6, 0, 1, 0, 1, 0, 0, 0, 0, 0] # expected Hessian - expected_symH = [2, 0, 0, 1, 0, 0, 0, 0, 6, 0, 0, 1, 0, 0, 0, 0, 0, 1] # expected symmetric Hessian - expected_symH2 = [2, 0, 0, sqrt(2), 0, 0, 0, 0, 6, 0, 0, sqrt(2), 0, 0, 0, 0, 0, sqrt(2)] # expected symmetric Hessian - expected_curl3 = [1 / 2 - 6 / 4, 0, 1 / 4 - 1 / 4] - - ## define grid = a single non-refenrece triangle - xgrid = reference_domain(Tetrahedron3D) - xgrid[Coordinates][:, 2] = [2, 0, 0] - - ## define P2-Courant finite element space - FEType = H1P2{3, 3} - FES = FESpace{FEType}(xgrid) - show(devnull, FES) - - ## get midpoint quadrature rule for constants - qf = QuadratureRule{Float64, Tetrahedron3D}(0) + update_basis!(FEBE_grad, 1) + update_basis!(FEBE_div, 1) + + @assert size(FEBE_grad.cvals, 1) == length(expected_grad) + @assert size(FEBE_div.cvals, 1) == length(expected_div) + grad = zeros(Float64, 9) + div = zeros(Float64, 1) + eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) + eval_febe!(div, FEBE_div, Iu.entries[FES[CellDofs][:, 1]], 1) + error_grad = sqrt(sum((grad - expected_grad) .^ 2)) + error_div = sqrt(sum((div - expected_div) .^ 2)) + println("EG = Tetrahedron3D | $fetype | operator = Gradient | error = $error_grad") + println("EG = Tetrahedron3D | $fetype | operator = Divergence | error = $error_div") + + if fetype <: AbstractHdivFiniteElement + return maximum([error_curl3, error_grad, error_div]) + end - ## define FE basis Evaluator for Hessian - FEBE_curl3 = FEEvaluator(FES, Curl3D, qf) FEBE_L = FEEvaluator(FES, Laplacian, qf) FEBE_H = FEEvaluator(FES, Hessian, qf) FEBE_symH = FEEvaluator(FES, SymmetricHessian{1}, qf) @@ -229,18 +272,12 @@ function test_derivatives3D() update_basis!(FEBE_H, 1) update_basis!(FEBE_symH, 1) update_basis!(FEBE_symH2, 1) - update_basis!(FEBE_curl3, 1) - - ## interpolate quadratic testfunction - Iu = FEVector(FES) - interpolate!(Iu[1], testf) ## check if operator evals have the correct length @assert size(FEBE_L.cvals, 1) == length(expected_L) @assert size(FEBE_H.cvals, 1) == length(expected_H) @assert size(FEBE_symH.cvals, 1) == length(expected_symH) @assert size(FEBE_symH2.cvals, 1) == length(expected_symH2) - @assert size(FEBE_curl3.cvals, 1) == length(expected_curl3) ## eval 2nd order derivatives at only quadrature point 1 ## since function is quadratic this should be constant @@ -248,74 +285,19 @@ function test_derivatives3D() symH = zeros(Float64, 18) symH2 = zeros(Float64, 18) L = zeros(Float64, 3) - curl3 = zeros(Float64, 3) eval_febe!(L, FEBE_L, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(H, FEBE_H, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(symH, FEBE_symH, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(symH2, FEBE_symH2, Iu.entries[FES[CellDofs][:, 1]], 1) - eval_febe!(curl3, FEBE_curl3, Iu.entries[FES[CellDofs][:, 1]], 1) ## compute errors to expected values error_L = sqrt(sum((L - expected_L) .^ 2)) error_H = sqrt(sum((H - expected_H) .^ 2)) error_symH = sqrt(sum((symH - expected_symH) .^ 2)) error_symH2 = sqrt(sum((symH2 - expected_symH2) .^ 2)) - error_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2)) - println("EG = Tetrahedron3D | operator = Curl3 | error = $error_curl3") - println("EG = Tetrahedron3D | operator = Laplacian | error = $error_L") - println("EG = Tetrahedron3D | operator = Hessian | error = $error_H") - println("EG = Tetrahedron3D | operator = SymmetricHessian{1} | error = $error_symH") - println("EG = Tetrahedron3D | operator = SymmetricHessian{√2} | error = $error_symH2") - - return maximum([error_curl3, error_L, error_H, error_symH, error_symH2]) -end - - -function test_derivatives3D_hdiv() - ## define test function and expected operator evals - function testf(result, qpinfo) - x = qpinfo.x - result[1] = x[2] - result[2] = 3 * x[3] - return result[3] = x[2] - end - - ## expected values of operators in cell midpoint - expected_curl3 = [1 - 3, 0, -1] - - ## define grid = a single non-refenrece triangle - xgrid = reference_domain(Tetrahedron3D) - xgrid[Coordinates][:, 2] = [2, 0, 0] - - ## define P2-Courant finite element space - FES = FESpace{HDIVRT1{3}}(xgrid) - show(devnull, FES) - - ## get midpoint quadrature rule for constants - qf = QuadratureRule{Float64, Tetrahedron3D}(0) - - ## define FE basis Evaluator for Hessian - FEBE_curl3 = FEEvaluator(FES, Curl3D, qf) - - ## update on cell 1 - update_basis!(FEBE_curl3, 1) - - ## interpolate quadratic testfunction - Iu = FEVector(FES) - interpolate!(Iu[1], testf) - - ## check if operator evals have the correct length - @assert size(FEBE_curl3.cvals, 1) == length(expected_curl3) - - ## eval 2nd order derivatives at only quadrature point 1 - ## since function is quadratic this should be constant - curl3 = zeros(Float64, 3) - eval_febe!(curl3, FEBE_curl3, Iu.entries[FES[CellDofs][:, 1]], 1) - @info curl3 - - ## compute errors to expected values - error_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2)) - println("EG = Tetrahedron3D | operator = Curl3 | error = $error_curl3") - - return maximum([error_curl3]) + println("EG = Tetrahedron3D | $fetype | operator = Laplacian | error = $error_L") + println("EG = Tetrahedron3D | $fetype | operator = Hessian | error = $error_H") + println("EG = Tetrahedron3D | $fetype | operator = SymmetricHessian{1} | error = $error_symH") + println("EG = Tetrahedron3D | $fetype | operator = SymmetricHessian{√2} | error = $error_symH2") + return maximum([error_curl3, error_L, error_H, error_symH, error_symH2, error_grad, error_div]) end