From fb402d49a9c084e41c74c7e58e2c7f0804711c10 Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Thu, 21 May 2026 14:30:05 +0200 Subject: [PATCH 01/13] improved evaluation of Gradient of HDIV functions, added more operator tests --- src/feevaluator_hdiv.jl | 16 +++++-- test/test_operators.jl | 96 ++++++++++++++++++++++++++++++----------- 2 files changed, 83 insertions(+), 29 deletions(-) diff --git a/src/feevaluator_hdiv.jl b/src/feevaluator_hdiv.jl index 6419842..c4032f1 100644 --- a/src/feevaluator_hdiv.jl +++ b/src/feevaluator_hdiv.jl @@ -96,12 +96,20 @@ 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) + for 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 j in 1:size(L2GM, 2) + temp = 0 + for m in 1:size(L2GAinv, 2) + temp += L2GAinv[k, m] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] + end + for c in 1:size(L2GM, 1) + cvals[k + offsets[c], dof_i, i] += L2GM[c, j] * temp + end + end + 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 diff --git a/test/test_operators.jl b/test/test_operators.jl index a943e8e..eac1040 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -78,9 +78,10 @@ function test_derivatives2D() 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] - ## define grid = a single non-refenrece triangle - xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]') + ## 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 P2-Courant finite element space FES = FESpace{H1P2{2, 2}}(xgrid) @@ -95,6 +96,7 @@ function test_derivatives2D() FEBE_H = FEEvaluator(FES, Hessian, qf) FEBE_symH = FEEvaluator(FES, SymmetricHessian{1}, qf) FEBE_symH2 = FEEvaluator(FES, SymmetricHessian{sqrt(2)}, qf) + FEBE_grad = FEEvaluator(FES, Gradient, qf) ## update on cell 1 update_basis!(FEBE_L, 1) @@ -102,7 +104,7 @@ function test_derivatives2D() update_basis!(FEBE_symH, 1) update_basis!(FEBE_symH2, 1) update_basis!(FEBE_curl2, 1) - + update_basis!(FEBE_grad, 1) ## interpolate quadratic testfunction Iu = FEVector(FES) interpolate!(Iu[1], testf) @@ -113,6 +115,7 @@ function test_derivatives2D() @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) + @assert size(FEBE_grad.cvals, 1) == length(expected_grad) ## eval 2nd order derivatives at only quadrature point 1 ## since function is quadratic this should be constant @@ -121,11 +124,13 @@ function test_derivatives2D() symH2 = zeros(Float64, 6) L = zeros(Float64, 2) curl2 = zeros(Float64, 1) + grad = zeros(Float64, 4) 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) + eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) ## compute errors to expected values error_L = sqrt(sum((L - expected_L) .^ 2)) @@ -133,13 +138,15 @@ function test_derivatives2D() 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]) + error_grad = sqrt(sum((grad - expected_grad) .^ 2)) + println("EG = Triangle2D | H1 | operator = Curl2 | error = $error_curl2") + println("EG = Triangle2D | H1 | operator = Laplacian | error = $error_L") + println("EG = Triangle2D | H1 | operator = Hessian | error = $error_H") + println("EG = Triangle2D | H1 | operator = SymmetricHessian{1} | error = $error_symH") + println("EG = Triangle2D | H1 | operator = SymmetricHessian{√2} | error = $error_symH2") + println("EG = Triangle2D | H1 | operator = Gradient | error = $error_grad") + + return maximum([error_curl2, error_L, error_H, error_symH, error_symH2, error_grad]) end @@ -153,9 +160,11 @@ function test_derivatives2D_hdiv() ## expected values of operators in cell midpoint expected_curl2 = [1 / 3] + expected_grad = [0, 0, 1 / 3, 2] + expected_div = [2] ## define grid = a single non-refenrece triangle - xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]') + xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]') # midpoint = [0.0, 1 / 3] ## define P2-Courant finite element space FES = FESpace{HDIVBDM2{2}}(xgrid) @@ -166,9 +175,13 @@ function test_derivatives2D_hdiv() ## define FE basis Evaluator for Hessian FEBE_curl2 = FEEvaluator(FES, Curl2D, qf) + FEBE_grad = FEEvaluator(FES, Gradient, qf) + FEBE_div = FEEvaluator(FES, Divergence, qf) ## update on cell 1 update_basis!(FEBE_curl2, 1) + update_basis!(FEBE_grad, 1) + update_basis!(FEBE_div, 1) ## interpolate quadratic testfunction Iu = FEVector(FES) @@ -176,16 +189,26 @@ function test_derivatives2D_hdiv() ## check if operator evals have the correct length @assert size(FEBE_curl2.cvals, 1) == length(expected_curl2) + @assert size(FEBE_grad.cvals, 1) == length(expected_grad) + @assert size(FEBE_div.cvals, 1) == length(expected_div) # compute derivatives curl2 = zeros(Float64, 1) + grad = zeros(Float64, 4) + div = zeros(Float64, 1) eval_febe!(curl2, FEBE_curl2, Iu.entries[FES[CellDofs][:, 1]], 1) + eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) + eval_febe!(div, FEBE_div, Iu.entries[FES[CellDofs][:, 1]], 1) ## compute errors to expected values error_curl2 = sqrt(sum((curl2 - expected_curl2) .^ 2)) - println("EG = Triangle2D | operator = Curl2 | error = $error_curl2") + error_grad = sqrt(sum((grad - expected_grad) .^ 2)) + error_div = sqrt(sum((div - expected_div) .^ 2)) + println("EG = Triangle2D | HDIV | operator = Curl2 | error = $error_curl2") + println("EG = Triangle2D | HDIV | operator = Gradient | error = $error_grad") + println("EG = Triangle2D | HDIV | operator = Divergence | error = $error_div") - return maximum([error_curl2]) + return maximum([error_curl2, error_grad, error_div]) end @@ -204,10 +227,11 @@ function test_derivatives3D() 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 ## define grid = a single non-refenrece triangle xgrid = reference_domain(Tetrahedron3D) - xgrid[Coordinates][:, 2] = [2, 0, 0] + xgrid[Coordinates][:, 2] = [2, 0, 0] # midpoint = [0.5, 0.25, 0.25] ## define P2-Courant finite element space FEType = H1P2{3, 3} @@ -223,6 +247,7 @@ function test_derivatives3D() FEBE_H = FEEvaluator(FES, Hessian, qf) FEBE_symH = FEEvaluator(FES, SymmetricHessian{1}, qf) FEBE_symH2 = FEEvaluator(FES, SymmetricHessian{sqrt(2)}, qf) + FEBE_grad = FEEvaluator(FES, Gradient, qf) ## update on cell 1 update_basis!(FEBE_L, 1) @@ -230,6 +255,7 @@ function test_derivatives3D() update_basis!(FEBE_symH, 1) update_basis!(FEBE_symH2, 1) update_basis!(FEBE_curl3, 1) + update_basis!(FEBE_grad, 1) ## interpolate quadratic testfunction Iu = FEVector(FES) @@ -241,6 +267,7 @@ function test_derivatives3D() @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) + @assert size(FEBE_grad.cvals, 1) == length(expected_grad) ## eval 2nd order derivatives at only quadrature point 1 ## since function is quadratic this should be constant @@ -249,11 +276,13 @@ function test_derivatives3D() symH2 = zeros(Float64, 18) L = zeros(Float64, 3) curl3 = zeros(Float64, 3) + grad = zeros(Float64, 9) 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) + eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) ## compute errors to expected values error_L = sqrt(sum((L - expected_L) .^ 2)) @@ -261,13 +290,15 @@ function test_derivatives3D() 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]) + error_grad = sqrt(sum((grad - expected_grad) .^ 2)) + println("EG = Tetrahedron3D | H1 | operator = Curl3 | error = $error_curl3") + println("EG = Tetrahedron3D | H1 |operator = Laplacian | error = $error_L") + println("EG = Tetrahedron3D | H1 | operator = Hessian | error = $error_H") + println("EG = Tetrahedron3D | H1 | operator = SymmetricHessian{1} | error = $error_symH") + println("EG = Tetrahedron3D | H1 | operator = SymmetricHessian{√2} | error = $error_symH2") + println("EG = Tetrahedron3D | H1 | operator = Gradient | error = $error_grad") + + return maximum([error_curl3, error_L, error_H, error_symH, error_symH2, error_grad]) end @@ -282,10 +313,12 @@ function test_derivatives3D_hdiv() ## expected values of operators in cell midpoint expected_curl3 = [1 - 3, 0, -1] + expected_grad = [0, 1, 0, 0, 0, 3, 0, 1, 0] # expected Gradient + expected_div = [0] ## define grid = a single non-refenrece triangle xgrid = reference_domain(Tetrahedron3D) - xgrid[Coordinates][:, 2] = [2, 0, 0] + xgrid[Coordinates][:, 2] = [2, 0, 0] # midpoint = [0.5, 0.25, 0.25] ## define P2-Courant finite element space FES = FESpace{HDIVRT1{3}}(xgrid) @@ -296,9 +329,13 @@ function test_derivatives3D_hdiv() ## define FE basis Evaluator for Hessian FEBE_curl3 = FEEvaluator(FES, Curl3D, qf) + FEBE_grad = FEEvaluator(FES, Gradient, qf) + FEBE_div = FEEvaluator(FES, Divergence, qf) ## update on cell 1 update_basis!(FEBE_curl3, 1) + update_basis!(FEBE_grad, 1) + update_basis!(FEBE_div, 1) ## interpolate quadratic testfunction Iu = FEVector(FES) @@ -306,16 +343,25 @@ function test_derivatives3D_hdiv() ## check if operator evals have the correct length @assert size(FEBE_curl3.cvals, 1) == length(expected_curl3) + @assert size(FEBE_grad.cvals, 1) == length(expected_grad) + @assert size(FEBE_div.cvals, 1) == length(expected_div) ## eval 2nd order derivatives at only quadrature point 1 ## since function is quadratic this should be constant curl3 = zeros(Float64, 3) + grad = zeros(Float64, 9) + div = zeros(Float64, 1) eval_febe!(curl3, FEBE_curl3, Iu.entries[FES[CellDofs][:, 1]], 1) - @info curl3 + eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) + eval_febe!(div, FEBE_div, Iu.entries[FES[CellDofs][:, 1]], 1) ## compute errors to expected values error_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2)) - println("EG = Tetrahedron3D | operator = Curl3 | error = $error_curl3") + error_grad = sqrt(sum((grad - expected_grad) .^ 2)) + error_div = sqrt(sum((div - expected_div) .^ 2)) + println("EG = Tetrahedron3D | HIDV | operator = Curl3 | error = $error_curl3") + println("EG = Tetrahedron3D | HIDV | operator = Gradient | error = $error_grad") + println("EG = Tetrahedron3D | HIDV | operator = Divergence | error = $error_div") - return maximum([error_curl3]) + return maximum([error_curl3, error_grad, error_div]) end From 14c5a4fb71e91587da4ffb0eaccd49753cbdaf57 Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Thu, 21 May 2026 14:48:54 +0200 Subject: [PATCH 02/13] improved curl evaluations for Hdiv --- src/feevaluator_h1.jl | 2 +- src/feevaluator_hdiv.jl | 32 ++++++++++++++++++++------------ 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/src/feevaluator_h1.jl b/src/feevaluator_h1.jl index ae6c7be..3435dbe 100644 --- a/src/feevaluator_h1.jl +++ b/src/feevaluator_h1.jl @@ -167,7 +167,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 diff --git a/src/feevaluator_hdiv.jl b/src/feevaluator_hdiv.jl index c4032f1..d0f018c 100644 --- a/src/feevaluator_hdiv.jl +++ b/src/feevaluator_hdiv.jl @@ -127,10 +127,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] + B = L2GAinv[1, m] * L2GM[2, j] + 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] / det + cvals[1, dof_i, i] += B * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det end end return nothing @@ -148,14 +150,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 From 4eb8056a41c0abefde110d97cfc81a3bb5f58779 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 13:35:16 +0200 Subject: [PATCH 03/13] further improvements of feevalutor loops (H1): dot is now used everywhere for application of inverse trafo to compute derivatives; removed some type annotations from get_coefficient closure functions; item loops and citem properties of structs now use Ti type --- src/fedefs/h1v_br.jl | 4 ++-- src/fedefs/hdiv_bdm1.jl | 4 ++-- src/fedefs/hdiv_bdm2.jl | 2 +- src/fedefs/hdiv_rtk.jl | 2 +- src/feevaluator.jl | 4 ++-- src/feevaluator_h1.jl | 38 ++++++++++++---------------------- src/feevaluator_hdiv.jl | 17 +++++++-------- src/interpolators.jl | 4 ++-- src/lazy_interpolate.jl | 8 +++---- src/point_evaluator.jl | 23 ++++++++++---------- src/quadrature.jl | 8 +++---- src/reconstructionoperators.jl | 2 +- 12 files changed, 51 insertions(+), 65 deletions(-) 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 3435dbe..f56c7c8 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] = dot(view(L2GAinv, k, :), view(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] = dot(view(L2GAinv, k, :), view(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 = dot(view(L2GAinv, k, :), view(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] += dot(view(L2GAinv, k, :), view(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] += dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[k], :, i)) * coefficients[k, dof_i] end return nothing end @@ -363,10 +353,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 = dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[c], :, i)) + cvals[1, dof_i, i] += temp * tangent[c] end end end diff --git a/src/feevaluator_hdiv.jl b/src/feevaluator_hdiv.jl index d0f018c..4c6a41a 100644 --- a/src/feevaluator_hdiv.jl +++ b/src/feevaluator_hdiv.jl @@ -97,16 +97,15 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad det = FEBE.L2G.det # 1 alloc for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) for k in 1:size(L2GAinv, 1) - # compute duc/dxk for j in 1:size(L2GM, 2) - temp = 0 - for m in 1:size(L2GAinv, 2) - temp += L2GAinv[k, m] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] - end + # apply inverse transform to d/dx_k of j-th component of i-th reference basis function + temp = dot(view(L2GAinv, k, :), view(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 @@ -128,11 +127,11 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl fill!(cvals, 0) det = FEBE.L2G.det # 1 alloc for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2) - A = L2GAinv[2, m] * L2GM[1, j] - B = L2GAinv[1, m] * L2GM[2, j] + 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] / det - cvals[1, dof_i, i] += B * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det + 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 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/point_evaluator.jl b/src/point_evaluator.jl index 2aa5743..28f9763 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 - ) + item::Ti, # 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, - ) + item::Ti, + ) where {Tv, Ti} ## find cell geometry id j = PE.eval_selector(item) 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) From 27e724f613142be59e3cd79c985771603dc03ba9 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 14:06:29 +0200 Subject: [PATCH 04/13] fixed point evaluator --- src/point_evaluator.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/point_evaluator.jl b/src/point_evaluator.jl index 28f9763..852ecfa 100644 --- a/src/point_evaluator.jl +++ b/src/point_evaluator.jl @@ -173,7 +173,7 @@ function initialize!(O::PointEvaluator{T, Ti, TCoeff, UT}, sol; time = 0, kwargs BE_args::Array{<:FEEvaluator, 1}, L2G::L2GTransformer{Tv, Ti}, xref, - item::Ti, # cell used to evaluate local coordinates + item, # cell used to evaluate local coordinates ) where {Tv, Ti} for id in 1:nargs @@ -247,7 +247,7 @@ function evaluate_bary!( result, PE::PointEvaluator{Tv, Ti}, xref, - item::Ti, + item, ) where {Tv, Ti} ## find cell geometry id From fedfc79cb18b05e4dc72744c6eca952312e9fa23 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 14:41:32 +0200 Subject: [PATCH 05/13] merged operator tests for H1 and Hdiv --- test/test_operators.jl | 340 +++++++++++++++-------------------------- 1 file changed, 125 insertions(+), 215 deletions(-) diff --git a/test/test_operators.jl b/test/test_operators.jl index eac1040..3dfb931 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -4,13 +4,15 @@ function run_operator_tests() println("============================") println("Testing Operator Evaluations") println("============================") - error = test_derivatives2D() + error = test_derivatives2D(H1P2{2, 2}) @test error < 1.0e-14 - error = test_derivatives2D_hdiv() + error = test_derivatives2D(HDIVBDM2{2}) @test error < 1.0e-14 - error = test_derivatives3D() + #error = test_derivatives2D_hdiv() + #@test error < 1.0e-14 + error = test_derivatives3D(H1P2{3, 3}, 2) @test error < 1.0e-14 - error = test_derivatives3D_hdiv() + error = test_derivatives3D(HDIVRT1{3}, 1) @test error < 1.0e-14 test_reconstructions() end @@ -64,7 +66,7 @@ function test_reconstructions() return end -function test_derivatives2D() +function test_derivatives2D(fetype) ## define test function and expected operator evals function testf(result, qpinfo) x = qpinfo.x @@ -72,6 +74,9 @@ function test_derivatives2D() return result[2] = 3 * x[2]^2 + x[1] * x[2] end + ## 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] + ## 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 @@ -79,97 +84,16 @@ function test_derivatives2D() expected_symH2 = [2, 0, 0, 0, 6, sqrt(2)] # expected symmetric Hessian expected_curl2 = [1 / 3] expected_grad = [0, 0, 1 / 3, 2] - - ## 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] + expected_div = [2] ## define P2-Courant finite element space - FES = FESpace{H1P2{2, 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) - 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) - FEBE_grad = FEEvaluator(FES, Gradient, 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) - update_basis!(FEBE_grad, 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) - @assert size(FEBE_grad.cvals, 1) == length(expected_grad) - - ## 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) - grad = zeros(Float64, 4) - 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) - eval_febe!(grad, FEBE_grad, 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)) - error_grad = sqrt(sum((grad - expected_grad) .^ 2)) - println("EG = Triangle2D | H1 | operator = Curl2 | error = $error_curl2") - println("EG = Triangle2D | H1 | operator = Laplacian | error = $error_L") - println("EG = Triangle2D | H1 | operator = Hessian | error = $error_H") - println("EG = Triangle2D | H1 | operator = SymmetricHessian{1} | error = $error_symH") - println("EG = Triangle2D | H1 | operator = SymmetricHessian{√2} | error = $error_symH2") - println("EG = Triangle2D | H1 | operator = Gradient | error = $error_grad") - - return maximum([error_curl2, error_L, error_H, error_symH, error_symH2, error_grad]) -end - - -function test_derivatives2D_hdiv() - ## 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] - end - - ## expected values of operators in cell midpoint - expected_curl2 = [1 / 3] - expected_grad = [0, 0, 1 / 3, 2] - expected_div = [2] - - ## define grid = a single non-refenrece triangle - xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]') # midpoint = [0.0, 1 / 3] - - ## define P2-Courant finite element space - FES = FESpace{HDIVBDM2{2}}(xgrid) - show(devnull, FES) - ## get midpoint quadrature rule for constants qf = QuadratureRule{Float64, Triangle2D}(0) @@ -177,191 +101,177 @@ function test_derivatives2D_hdiv() FEBE_curl2 = FEEvaluator(FES, Curl2D, qf) FEBE_grad = FEEvaluator(FES, Gradient, qf) FEBE_div = FEEvaluator(FES, Divergence, qf) - - ## update on cell 1 update_basis!(FEBE_curl2, 1) update_basis!(FEBE_grad, 1) update_basis!(FEBE_div, 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) @assert size(FEBE_grad.cvals, 1) == length(expected_grad) - @assert size(FEBE_div.cvals, 1) == length(expected_div) - - # compute derivatives + # evaluate at quadrature points = cell midpoint curl2 = zeros(Float64, 1) grad = zeros(Float64, 4) div = zeros(Float64, 1) eval_febe!(curl2, FEBE_curl2, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(div, FEBE_div, Iu.entries[FES[CellDofs][:, 1]], 1) - ## compute errors to expected values error_curl2 = sqrt(sum((curl2 - expected_curl2) .^ 2)) error_grad = sqrt(sum((grad - expected_grad) .^ 2)) error_div = sqrt(sum((div - expected_div) .^ 2)) - println("EG = Triangle2D | HDIV | operator = Curl2 | error = $error_curl2") - println("EG = Triangle2D | HDIV | operator = Gradient | error = $error_grad") - println("EG = Triangle2D | HDIV | operator = Divergence | error = $error_div") + println("EG = Triangle2D | $fetype | operator = Curl2 | error = $error_curl2") + println("EG = Triangle2D | $fetype | operator = Gradient | error = $error_grad") + println("EG = Triangle2D | $fetype | operator = Divergence | error = $error_div") + + if fetype <: AbstractH1FiniteElement + # 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_basis!(FEBE_L, 1) + update_basis!(FEBE_H, 1) + update_basis!(FEBE_symH, 1) + update_basis!(FEBE_symH2, 1) + @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) + H = zeros(Float64, 8) + symH = zeros(Float64, 6) + symH2 = zeros(Float64, 6) + L = zeros(Float64, 2) + 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) + 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)) + 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_curl2, error_L, error_H, error_symH, error_symH2, error_grad, error_div]) + else + return maximum([error_curl2, error_grad]) + end - return maximum([error_curl2, error_grad, error_div]) -end +end -function test_derivatives3D() +function test_derivatives3D(fetype, order) ## 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 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_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 - ## define grid = a single non-refenrece triangle xgrid = reference_domain(Tetrahedron3D) xgrid[Coordinates][:, 2] = [2, 0, 0] # midpoint = [0.5, 0.25, 0.25] ## define P2-Courant finite element space - FEType = H1P2{3, 3} - FES = FESpace{FEType}(xgrid) + FES = FESpace{fetype}(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) - 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) - FEBE_grad = FEEvaluator(FES, Gradient, 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_curl3, 1) - update_basis!(FEBE_grad, 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) - @assert size(FEBE_grad.cvals, 1) == length(expected_grad) - - ## eval 2nd order derivatives at only quadrature point 1 - ## since function is quadratic this should be constant - H = zeros(Float64, 27) - symH = zeros(Float64, 18) - symH2 = zeros(Float64, 18) - L = zeros(Float64, 3) - curl3 = zeros(Float64, 3) - grad = zeros(Float64, 9) - 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) - eval_febe!(grad, FEBE_grad, 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)) - error_grad = sqrt(sum((grad - expected_grad) .^ 2)) - println("EG = Tetrahedron3D | H1 | operator = Curl3 | error = $error_curl3") - println("EG = Tetrahedron3D | H1 |operator = Laplacian | error = $error_L") - println("EG = Tetrahedron3D | H1 | operator = Hessian | error = $error_H") - println("EG = Tetrahedron3D | H1 | operator = SymmetricHessian{1} | error = $error_symH") - println("EG = Tetrahedron3D | H1 | operator = SymmetricHessian{√2} | error = $error_symH2") - println("EG = Tetrahedron3D | H1 | operator = Gradient | error = $error_grad") - - return maximum([error_curl3, error_L, error_H, error_symH, error_symH2, error_grad]) -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] - expected_grad = [0, 1, 0, 0, 0, 3, 0, 1, 0] # expected Gradient - expected_div = [0] - - ## define grid = a single non-refenrece triangle - 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{HDIVRT1{3}}(xgrid) - show(devnull, FES) + 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 ## get midpoint quadrature rule for constants qf = QuadratureRule{Float64, Tetrahedron3D}(0) ## define FE basis Evaluator for Hessian - FEBE_curl3 = FEEvaluator(FES, Curl3D, qf) FEBE_grad = FEEvaluator(FES, Gradient, qf) FEBE_div = FEEvaluator(FES, Divergence, qf) + FEBE_curl3 = FEEvaluator(FES, Curl3D, qf) - ## update on cell 1 update_basis!(FEBE_curl3, 1) update_basis!(FEBE_grad, 1) update_basis!(FEBE_div, 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) @assert size(FEBE_grad.cvals, 1) == length(expected_grad) @assert size(FEBE_div.cvals, 1) == length(expected_div) - - ## eval 2nd order derivatives at only quadrature point 1 - ## since function is quadratic this should be constant curl3 = zeros(Float64, 3) grad = zeros(Float64, 9) div = zeros(Float64, 1) eval_febe!(curl3, FEBE_curl3, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) eval_febe!(div, FEBE_div, Iu.entries[FES[CellDofs][:, 1]], 1) - - ## compute errors to expected values - error_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2)) error_grad = sqrt(sum((grad - expected_grad) .^ 2)) + error_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2)) error_div = sqrt(sum((div - expected_div) .^ 2)) - println("EG = Tetrahedron3D | HIDV | operator = Curl3 | error = $error_curl3") - println("EG = Tetrahedron3D | HIDV | operator = Gradient | error = $error_grad") - println("EG = Tetrahedron3D | HIDV | operator = Divergence | error = $error_div") + println("EG = Tetrahedron3D | $fetype | operator = Curl3 | error = $error_curl3") + println("EG = Tetrahedron3D | $fetype | operator = Gradient | error = $error_grad") + println("EG = Tetrahedron3D | $fetype | operator = Divergence | error = $error_div") + + if fetype <: AbstractH1FiniteElement + 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) + + ## 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) + + ## eval 2nd order derivatives at only quadrature point 1 + ## since function is quadratic this should be constant + H = zeros(Float64, 27) + symH = zeros(Float64, 18) + symH2 = zeros(Float64, 18) + L = 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) + + ## 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)) + 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]) + else + return maximum([error_curl3, error_grad, error_div]) + end - return maximum([error_curl3, error_grad, error_div]) end From fe2beaec9d09c0d9ad2e93ade4deacf2094e8d95 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 14:48:53 +0200 Subject: [PATCH 06/13] further improvements in curl operators for H1 --- src/feevaluator_h1.jl | 41 +++++++++++++---------------------------- 1 file changed, 13 insertions(+), 28 deletions(-) diff --git a/src/feevaluator_h1.jl b/src/feevaluator_h1.jl index f56c7c8..eb05508 100644 --- a/src/feevaluator_h1.jl +++ b/src/feevaluator_h1.jl @@ -278,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 + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] = -dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i], :, i)) # -du1/dy + cvals[2, dof_i, i] = dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i], :, i)) # du2/dx end return nothing end @@ -298,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 + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] = dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i] + offsets2[1], :, i)) # -du1/dy + cvals[1, dof_i, i] += dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i] + offsets2[2], :, i)) # du2/dx end return nothing end @@ -319,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 + for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) + cvals[1, dof_i, i] = dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i] + offsets2[3], :, i)) # du3/dx2 + cvals[1, dof_i, i] -= dot(view(L2GAinv, 3, :), view(refbasisderivvals, subset[dof_i] + offsets2[2], :, i)) # - du2/dx3 + cvals[2, dof_i, i] = dot(view(L2GAinv, 3, :), view(refbasisderivvals, subset[dof_i] + offsets2[1], :, i)) # du1/dx3 + cvals[2, dof_i, i] -= dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i] + offsets2[3], :, i)) # - du3/dx1 + cvals[3, dof_i, i] = dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i] + offsets2[2], :, i)) # du2/dx1 + cvals[3, dof_i, i] -= dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i] + offsets2[1], :, i)) # - du1/dx2 end return nothing end From e7ddf9894887f8f06cf34bf51a023f7d75ff5b54 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 15:10:20 +0200 Subject: [PATCH 07/13] curls of Hcurl transformation now also tested, fixed sign in Curl2 of Hdiv --- src/feevaluator_hcurl.jl | 4 +- test/test_operators.jl | 255 +++++++++++++++++++++------------------ 2 files changed, 140 insertions(+), 119 deletions(-) diff --git a/src/feevaluator_hcurl.jl b/src/feevaluator_hcurl.jl index f5a6806..4a2bbe6 100644 --- a/src/feevaluator_hcurl.jl +++ b/src/feevaluator_hcurl.jl @@ -43,8 +43,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/test/test_operators.jl b/test/test_operators.jl index 3dfb931..36b7890 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -4,16 +4,18 @@ function run_operator_tests() println("============================") println("Testing Operator Evaluations") println("============================") - error = test_derivatives2D(H1P2{2, 2}) + error = test_derivatives2D(H1P2{2, 2}, 2) @test error < 1.0e-14 - error = test_derivatives2D(HDIVBDM2{2}) + error = test_derivatives2D(HDIVBDM2{2}, 2) + @test error < 1.0e-14 + error = test_derivatives2D(HCURLN1{2}, 2) @test error < 1.0e-14 - #error = test_derivatives2D_hdiv() - #@test error < 1.0e-14 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 end @@ -66,25 +68,35 @@ function test_reconstructions() return end -function test_derivatives2D(fetype) +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 ## 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] ## 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] - expected_grad = [0, 0, 1 / 3, 2] - expected_div = [2] + if order == 2 + 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_curl2 = [-2] + end ## define P2-Courant finite element space FES = FESpace{fetype}(xgrid) @@ -97,64 +109,68 @@ function test_derivatives2D(fetype) ## get midpoint quadrature rule for constants qf = QuadratureRule{Float64, Triangle2D}(0) - ## define FE basis Evaluator for Hessian FEBE_curl2 = FEEvaluator(FES, Curl2D, qf) - FEBE_grad = FEEvaluator(FES, Gradient, qf) - FEBE_div = FEEvaluator(FES, Divergence, qf) update_basis!(FEBE_curl2, 1) - update_basis!(FEBE_grad, 1) - update_basis!(FEBE_div, 1) ## check if operator evals have the correct length @assert size(FEBE_curl2.cvals, 1) == length(expected_curl2) - @assert size(FEBE_grad.cvals, 1) == length(expected_grad) # evaluate at quadrature points = cell midpoint curl2 = zeros(Float64, 1) grad = zeros(Float64, 4) div = zeros(Float64, 1) eval_febe!(curl2, FEBE_curl2, Iu.entries[FES[CellDofs][:, 1]], 1) - eval_febe!(grad, FEBE_grad, Iu.entries[FES[CellDofs][:, 1]], 1) - eval_febe!(div, FEBE_div, Iu.entries[FES[CellDofs][:, 1]], 1) ## compute errors to expected values error_curl2 = sqrt(sum((curl2 - expected_curl2) .^ 2)) - error_grad = sqrt(sum((grad - expected_grad) .^ 2)) - error_div = sqrt(sum((div - expected_div) .^ 2)) println("EG = Triangle2D | $fetype | operator = Curl2 | error = $error_curl2") - println("EG = Triangle2D | $fetype | operator = Gradient | error = $error_grad") - println("EG = Triangle2D | $fetype | operator = Divergence | error = $error_div") - - if fetype <: AbstractH1FiniteElement - # 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_basis!(FEBE_L, 1) - update_basis!(FEBE_H, 1) - update_basis!(FEBE_symH, 1) - update_basis!(FEBE_symH2, 1) - @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) - H = zeros(Float64, 8) - symH = zeros(Float64, 6) - symH2 = zeros(Float64, 6) - L = zeros(Float64, 2) - 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) - 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)) - 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_curl2, error_L, error_H, error_symH, error_symH2, error_grad, error_div]) + if fetype <: AbstractHcurlFiniteElement + return maximum([error_curl2]) else - return maximum([error_curl2, error_grad]) + 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 <: AbstractH1FiniteElement + # 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_basis!(FEBE_L, 1) + update_basis!(FEBE_H, 1) + update_basis!(FEBE_symH, 1) + update_basis!(FEBE_symH2, 1) + @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) + H = zeros(Float64, 8) + symH = zeros(Float64, 6) + symH2 = zeros(Float64, 6) + L = zeros(Float64, 2) + 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) + 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)) + 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_curl2, error_L, error_H, error_symH, error_symH2, error_grad, error_div]) + else + return maximum([error_curl2, error_grad]) + end end @@ -206,72 +222,77 @@ function test_derivatives3D(fetype, order) ## get midpoint quadrature rule for constants qf = QuadratureRule{Float64, Tetrahedron3D}(0) - ## define FE basis Evaluator for Hessian - FEBE_grad = FEEvaluator(FES, Gradient, qf) - FEBE_div = FEEvaluator(FES, Divergence, qf) FEBE_curl3 = FEEvaluator(FES, Curl3D, qf) - update_basis!(FEBE_curl3, 1) - update_basis!(FEBE_grad, 1) - update_basis!(FEBE_div, 1) - @assert size(FEBE_curl3.cvals, 1) == length(expected_curl3) - @assert size(FEBE_grad.cvals, 1) == length(expected_grad) - @assert size(FEBE_div.cvals, 1) == length(expected_div) + + FEBE_grad = FEEvaluator(FES, Gradient, qf) + FEBE_div = FEEvaluator(FES, Divergence, qf) curl3 = zeros(Float64, 3) - grad = zeros(Float64, 9) - div = zeros(Float64, 1) eval_febe!(curl3, FEBE_curl3, Iu.entries[FES[CellDofs][:, 1]], 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_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2)) - error_div = sqrt(sum((div - expected_div) .^ 2)) println("EG = Tetrahedron3D | $fetype | operator = Curl3 | error = $error_curl3") - println("EG = Tetrahedron3D | $fetype | operator = Gradient | error = $error_grad") - println("EG = Tetrahedron3D | $fetype | operator = Divergence | error = $error_div") - - if fetype <: AbstractH1FiniteElement - 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) - - ## 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) - - ## eval 2nd order derivatives at only quadrature point 1 - ## since function is quadratic this should be constant - H = zeros(Float64, 27) - symH = zeros(Float64, 18) - symH2 = zeros(Float64, 18) - L = 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) - - ## 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)) - 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]) + + if fetype <: AbstractHcurlFiniteElement + return maximum([error_curl3]) else - return maximum([error_curl3, error_grad, error_div]) + + 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 <: AbstractH1FiniteElement + 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) + + ## 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) + + ## eval 2nd order derivatives at only quadrature point 1 + ## since function is quadratic this should be constant + H = zeros(Float64, 27) + symH = zeros(Float64, 18) + symH2 = zeros(Float64, 18) + L = 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) + + ## 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)) + 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]) + else + return maximum([error_curl3, error_grad, error_div]) + end end end From d5322b541aaa6e1b51e1520c561ae02c12d07bca Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 15:36:20 +0200 Subject: [PATCH 08/13] version bump + changelog --- CHANGELOG.md | 3 +++ Project.toml | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d200572..8ca9a27 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,8 @@ # CHANGES +## v1.7.0 May 22, 2026 +- improved operator evaluations, more operator tests + ## 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] From 100a127ee6255f4325d82df0a2a65963b093c153 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 16:49:17 +0200 Subject: [PATCH 09/13] improved Identity evaluations for Hdiv, Hcurl, added 2D Identity test --- src/feevaluator_hcurl.jl | 6 +----- src/feevaluator_hdiv.jl | 12 ++---------- test/test_operators.jl | 21 +++++++++++++++------ 3 files changed, 18 insertions(+), 21 deletions(-) diff --git a/src/feevaluator_hcurl.jl b/src/feevaluator_hcurl.jl index 4a2bbe6..f14490b 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] = dot(view(L2GAinv, k, :), view(refbasisvals[i], subset[dof_i], :)) * coefficients[k, dof_i] end return nothing end diff --git a/src/feevaluator_hdiv.jl b/src/feevaluator_hdiv.jl index 4c6a41a..310f2d0 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] = dot(view(L2GM, k, :), view(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] = dot(view(L2GM, c, :), view(refbasisvals[i], subset[dof_i], :)) * coefficients[c, dof_i] / det end return nothing end diff --git a/test/test_operators.jl b/test/test_operators.jl index 36b7890..191e05e 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -8,7 +8,7 @@ function run_operator_tests() @test error < 1.0e-14 error = test_derivatives2D(HDIVBDM2{2}, 2) @test error < 1.0e-14 - error = test_derivatives2D(HCURLN1{2}, 2) + error = test_derivatives2D(HCURLN1{2}, 1) @test error < 1.0e-14 error = test_derivatives3D(H1P2{3, 3}, 2) @test error < 1.0e-14 @@ -87,6 +87,7 @@ function test_derivatives2D(fetype, order) ## 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 @@ -95,6 +96,7 @@ function test_derivatives2D(fetype, order) expected_grad = [0, 0, 1 / 3, 2] expected_div = [2] elseif order == 1 + expected_id = [4 / 3, 1] expected_curl2 = [-2] end @@ -109,21 +111,28 @@ function test_derivatives2D(fetype, order) ## get midpoint quadrature rule for constants qf = QuadratureRule{Float64, Triangle2D}(0) + 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) - grad = zeros(Float64, 4) - div = 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_curl2]) + return maximum([error_id, error_curl2]) else + 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) @@ -167,9 +176,9 @@ function test_derivatives2D(fetype, order) 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_curl2, error_L, error_H, error_symH, error_symH2, error_grad, error_div]) + return maximum([error_id, error_curl2, error_L, error_H, error_symH, error_symH2, error_grad, error_div]) else - return maximum([error_curl2, error_grad]) + return maximum([error_id, error_curl2, error_grad]) end end From d9fbe6e04478c93c65afa406cf19f2590ad92082 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 17:30:21 +0200 Subject: [PATCH 10/13] use tags in Examples, fixed Example210: errors look ok now, scalarplots for VectorBlock now has a default title --- examples/Example200_LowLevelPoisson.jl | 14 +++---- .../Example205_LowLevelSpaceTimePoisson.jl | 14 +++---- examples/Example210_LowLevelNavierStokes.jl | 37 +++++++++++-------- src/plots.jl | 14 +++---- src/qpinfos.jl | 2 +- 5 files changed, 43 insertions(+), 38 deletions(-) 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..bfb57bd 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,8 +208,8 @@ 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 - update_trafo!(L2G_time, time_cell) + for time_cell::Ti in 1:ncells_time + update_trafo!(L2G_time, time_cell, Val(true)) for jT in 1:ndofs4cell_time, kT in 1:ndofs4cell_time dofTj = celldofs_time[jT, time_cell] dofTk = celldofs_time[kT, time_cell] @@ -233,12 +233,12 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f, fill!(Aloc, 0) ## assemble right-hand side - update_trafo!(L2G_space, cell) + update_trafo!(L2G_space, cell, Val(true)) 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 - update_trafo!(L2G_time, time_cell) + for time_cell::Ti in 1:ncells_time + update_trafo!(L2G_time, time_cell, Val(true)) for qpT in 1:nweights_time ## evaluate time coordinate eval_trafo!(t, L2G_time, xref_time[qpT]) 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/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/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. From 8965633ed1c6527b5f1dd2b546bff383549a622b Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 18:01:16 +0200 Subject: [PATCH 11/13] views macro and improvements suggested by @pjaap --- src/feevaluator_h1.jl | 38 +++---- src/feevaluator_hcurl.jl | 2 +- src/feevaluator_hdiv.jl | 6 +- test/test_operators.jl | 212 +++++++++++++++++++-------------------- 4 files changed, 127 insertions(+), 131 deletions(-) diff --git a/src/feevaluator_h1.jl b/src/feevaluator_h1.jl index eb05508..57640fb 100644 --- a/src/feevaluator_h1.jl +++ b/src/feevaluator_h1.jl @@ -63,7 +63,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad refbasisderivvals = FEBE.refbasisderivvals 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] = dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[c], :, i)) + cvals[k + offsets[c], dof_i, i] = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, i]) end return nothing end @@ -78,7 +78,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals 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) - cvals[k + offsets[c], dof_i, i] = dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[c], :, 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 @@ -95,7 +95,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Symm 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) # compute duc/dxk and put it into the right spot in the Voigt vector - temp = dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[c], :, i)) + temp = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, i]) if k != c cvals[compression[k + offsets[c]], dof_i, i] += offdiagval * temp else @@ -114,7 +114,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Dive 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) - cvals[1, dof_i, i] += dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[k], :, i)) + cvals[1, dof_i, i] += @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[k], :, i]) end return nothing end @@ -129,7 +129,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Dive 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) - cvals[1, dof_i, i] += dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[k], :, i)) * coefficients[k, dof_i] + cvals[1, dof_i, i] += @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[k], :, i]) * coefficients[k, dof_i] end return nothing end @@ -278,9 +278,9 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - cvals[1, dof_i, i] = -dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i], :, i)) # -du1/dy - cvals[2, dof_i, i] = dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i], :, i)) # du2/dx + @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 @@ -293,9 +293,9 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - cvals[1, dof_i, i] = dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i] + offsets2[1], :, i)) # -du1/dy - cvals[1, dof_i, i] += dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i] + offsets2[2], :, i)) # du2/dx + @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 @@ -309,13 +309,13 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl offsets = FEBE.offsets offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals - for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - cvals[1, dof_i, i] = dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i] + offsets2[3], :, i)) # du3/dx2 - cvals[1, dof_i, i] -= dot(view(L2GAinv, 3, :), view(refbasisderivvals, subset[dof_i] + offsets2[2], :, i)) # - du2/dx3 - cvals[2, dof_i, i] = dot(view(L2GAinv, 3, :), view(refbasisderivvals, subset[dof_i] + offsets2[1], :, i)) # du1/dx3 - cvals[2, dof_i, i] -= dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i] + offsets2[3], :, i)) # - du3/dx1 - cvals[3, dof_i, i] = dot(view(L2GAinv, 1, :), view(refbasisderivvals, subset[dof_i] + offsets2[2], :, i)) # du2/dx1 - cvals[3, dof_i, i] -= dot(view(L2GAinv, 2, :), view(refbasisderivvals, subset[dof_i] + offsets2[1], :, i)) # - du1/dx2 + @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 @@ -338,7 +338,7 @@ 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) - temp = dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[c], :, i)) + temp = @views dot(L2GAinv[k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, i]) cvals[1, dof_i, i] += temp * tangent[c] end end diff --git a/src/feevaluator_hcurl.jl b/src/feevaluator_hcurl.jl index f14490b..67ea42b 100644 --- a/src/feevaluator_hcurl.jl +++ b/src/feevaluator_hcurl.jl @@ -6,7 +6,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden refbasisvals = FEBE.refbasisvals cvals = FEBE.cvals 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] = dot(view(L2GAinv, k, :), view(refbasisvals[i], subset[dof_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 diff --git a/src/feevaluator_hdiv.jl b/src/feevaluator_hdiv.jl index 310f2d0..62ac927 100644 --- a/src/feevaluator_hdiv.jl +++ b/src/feevaluator_hdiv.jl @@ -8,7 +8,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden refbasisvals = FEBE.refbasisvals for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) for k in 1:size(L2GM, 1) - cvals[k, dof_i, i] = dot(view(L2GM, k, :), view(refbasisvals[i], subset[dof_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 @@ -24,7 +24,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Iden cvals = FEBE.cvals refbasisvals = FEBE.refbasisvals for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2) - cvals[1, dof_i, i] = dot(view(L2GM, c, :), view(refbasisvals[i], subset[dof_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 @@ -91,7 +91,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad 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 = dot(view(L2GAinv, k, :), view(refbasisderivvals, subset[dof_i] + offsets2[j], :, i)) + 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 diff --git a/test/test_operators.jl b/test/test_operators.jl index 191e05e..48c2c71 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -130,59 +130,57 @@ function test_derivatives2D(fetype, order) println("EG = Triangle2D | $fetype | operator = Curl2 | error = $error_curl2") if fetype <: AbstractHcurlFiniteElement return maximum([error_id, error_curl2]) - else - 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 <: AbstractH1FiniteElement - # 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_basis!(FEBE_L, 1) - update_basis!(FEBE_H, 1) - update_basis!(FEBE_symH, 1) - update_basis!(FEBE_symH2, 1) - @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) - H = zeros(Float64, 8) - symH = zeros(Float64, 6) - symH2 = zeros(Float64, 6) - L = zeros(Float64, 2) - 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) - 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)) - 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]) - else - return maximum([error_id, error_curl2, error_grad]) - end 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_basis!(FEBE_L, 1) + update_basis!(FEBE_H, 1) + update_basis!(FEBE_symH, 1) + update_basis!(FEBE_symH2, 1) + @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) + H = zeros(Float64, 8) + symH = zeros(Float64, 6) + symH2 = zeros(Float64, 6) + L = zeros(Float64, 2) + 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) + 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)) + 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_derivatives3D(fetype, order) @@ -244,64 +242,62 @@ function test_derivatives3D(fetype, order) if fetype <: AbstractHcurlFiniteElement return maximum([error_curl3]) - else - - 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 <: AbstractH1FiniteElement - 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) - - ## 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) - - ## eval 2nd order derivatives at only quadrature point 1 - ## since function is quadratic this should be constant - H = zeros(Float64, 27) - symH = zeros(Float64, 18) - symH2 = zeros(Float64, 18) - L = 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) - - ## 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)) - 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]) - else - return maximum([error_curl3, error_grad, error_div]) - end end + 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 + + 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) + + ## 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) + + ## eval 2nd order derivatives at only quadrature point 1 + ## since function is quadratic this should be constant + H = zeros(Float64, 27) + symH = zeros(Float64, 18) + symH2 = zeros(Float64, 18) + L = 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) + + ## 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)) + 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 From 83e32ed3c1bd83e3c666b5900336b94991a227c8 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Fri, 22 May 2026 19:02:11 +0200 Subject: [PATCH 12/13] fixed bug introduced in last commit, updated changelog --- CHANGELOG.md | 4 +++- src/feevaluator_h1.jl | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8ca9a27..97cc43e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,9 @@ # CHANGES -## v1.7.0 May 22, 2026 +## 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/src/feevaluator_h1.jl b/src/feevaluator_h1.jl index 57640fb..963dd1d 100644 --- a/src/feevaluator_h1.jl +++ b/src/feevaluator_h1.jl @@ -78,7 +78,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad offsets2 = FEBE.offsets2 refbasisderivvals = FEBE.refbasisderivvals 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) - cvals[k + offsets[c], dof_i, i] = @views dot([L2GAinv, k, :], refbasisderivvals[subset[dof_i] + offsets2[c], :, 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 From cd04ddda3703f9afd696ad5692e9aa80b7a2a86d Mon Sep 17 00:00:00 2001 From: chmerdon Date: Sat, 23 May 2026 17:25:32 +0200 Subject: [PATCH 13/13] fixed Example205 --- examples/Example205_LowLevelSpaceTimePoisson.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/Example205_LowLevelSpaceTimePoisson.jl b/examples/Example205_LowLevelSpaceTimePoisson.jl index bfb57bd..a6de045 100644 --- a/examples/Example205_LowLevelSpaceTimePoisson.jl +++ b/examples/Example205_LowLevelSpaceTimePoisson.jl @@ -209,7 +209,7 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f, ## add local matrix to global matrix for time_cell::Ti in 1:ncells_time - update_trafo!(L2G_time, time_cell, Val(true)) + update_trafo!(L2G_time, time_cell) for jT in 1:ndofs4cell_time, kT in 1:ndofs4cell_time dofTj = celldofs_time[jT, time_cell] dofTk = celldofs_time[kT, time_cell] @@ -233,12 +233,12 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f, fill!(Aloc, 0) ## assemble right-hand side - update_trafo!(L2G_space, cell, Val(true)) + update_trafo!(L2G_space, cell) for qp in 1:nweights_space ## evaluate coordinates of quadrature point in space eval_trafo!(x, L2G_space, xref_space[qp]) for time_cell::Ti in 1:ncells_time - update_trafo!(L2G_time, time_cell, Val(true)) + update_trafo!(L2G_time, time_cell) for qpT in 1:nweights_time ## evaluate time coordinate eval_trafo!(t, L2G_time, xref_time[qpT])