From a53680ff4779e2ed5751c8a922ea9856c5a85e0b Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Thu, 15 May 2025 12:07:50 +0200 Subject: [PATCH] Add Runic Code Quality check and badge --- .pre-commit-config.yaml | 4 + README.md | 1 + docs/make.jl | 68 +- examples/Example200_LowLevelPoisson.jl | 316 ++++----- .../Example205_LowLevelSpaceTimePoisson.jl | 426 ++++++------ examples/Example210_LowLevelNavierStokes.jl | 650 +++++++++--------- examples/Example280_BasisPlotter.jl | 36 +- examples/Example281_DiscontinuousPlot.jl | 38 +- .../Example290_InterpolationBetweenMeshes.jl | 60 +- pluto-examples/LowLevelNavierStokes.jl | 554 +++++++-------- pluto-examples/LowLevelPoisson.jl | 217 +++--- src/fedefs/h1_bubble.jl | 2 +- src/fedefs/h1_p2b.jl | 2 +- src/fedefs/h1v_p1teb.jl | 2 +- src/fedefs/hcurl_n0.jl | 2 +- src/fedefs/hdiv_bdm1.jl | 6 +- src/fedefs/hdiv_bdm2.jl | 10 +- src/fedefs/hdiv_rt1.jl | 6 +- src/fedefs/hdiv_rtk.jl | 5 +- src/fematrix.jl | 15 +- src/fevector.jl | 8 +- src/finiteelements.jl | 2 +- src/interpolators.jl | 98 +-- test/test_interpolators.jl | 6 +- test/test_segmentintegrator.jl | 2 +- 25 files changed, 1276 insertions(+), 1260 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b28f63a0..14f576c7 100755 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -19,3 +19,7 @@ repos: rev: v2.4.1 hooks: - id: codespell +- repo: https://github.com/fredrikekre/runic-pre-commit + rev: v1.0.0 + hooks: + - id: runic diff --git a/README.md b/README.md index 29790970..c9c23fb8 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,7 @@ [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://wias-pdelib.github.io/ExtendableFEMBase.jl/stable/index.html) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://wias-pdelib.github.io/ExtendableFEMBase.jl/dev/index.html) [![DOI](https://zenodo.org/badge/667751152.svg)](https://zenodo.org/doi/10.5281/zenodo.10563410) +[![code style: runic](https://img.shields.io/badge/code_style-%E1%9A%B1%E1%9A%A2%E1%9A%BE%E1> # ExtendableFEMBase diff --git a/docs/make.jl b/docs/make.jl index 3cbb15bb..a41f7638 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,42 +7,42 @@ using CairoMakie function make_all(; with_examples::Bool = true) - module_examples = [] - pluto_examples = [] + module_examples = [] + pluto_examples = [] if with_examples - - DocMeta.setdocmeta!(ExampleJuggler, :DocTestSetup, :(using ExampleJuggler); recursive = true) - - example_dir = joinpath(@__DIR__, "..", "examples") - pluto_example_dir = joinpath(@__DIR__, "..", "pluto-examples") - - modules = [ - "Example200_LowLevelPoisson.jl", - "Example205_LowLevelSpaceTimePoisson.jl", - "Example210_LowLevelNavierStokes.jl", - "Example280_BasisPlotter.jl", - "Example281_DiscontinuousPlot.jl", - "Example290_InterpolationBetweenMeshes.jl", - ] - + + DocMeta.setdocmeta!(ExampleJuggler, :DocTestSetup, :(using ExampleJuggler); recursive = true) + + example_dir = joinpath(@__DIR__, "..", "examples") + pluto_example_dir = joinpath(@__DIR__, "..", "pluto-examples") + + modules = [ + "Example200_LowLevelPoisson.jl", + "Example205_LowLevelSpaceTimePoisson.jl", + "Example210_LowLevelNavierStokes.jl", + "Example280_BasisPlotter.jl", + "Example281_DiscontinuousPlot.jl", + "Example290_InterpolationBetweenMeshes.jl", + ] + notebooks = [ "Low level Poisson" => "LowLevelPoisson.jl" "Low level Navier-Stokes" => "LowLevelNavierStokes.jl" ] - cleanexamples() + cleanexamples() - module_examples = @docmodules(example_dir, modules, Plotter = CairoMakie) - pluto_examples = @docplutonotebooks(pluto_example_dir, notebooks, iframe=false) + module_examples = @docmodules(example_dir, modules, Plotter = CairoMakie) + pluto_examples = @docplutonotebooks(pluto_example_dir, notebooks, iframe = false) pushfirst!(module_examples, "Introduction" => "examples_intro.md") end makedocs( - modules=[ExtendableFEMBase], - sitename="ExtendableFEMBase.jl", - authors="Christian Merdon", + modules = [ExtendableFEMBase], + sitename = "ExtendableFEMBase.jl", + authors = "Christian Merdon", format = Documenter.HTML(repolink = "https://github.com/WIAS-PDELib/ExtendableFEMBase.jl", size_threshold = 250000, mathengine = MathJax3()), clean = false, checkdocs = :none, # :all or :exports currently causes UndefVarError @@ -53,26 +53,26 @@ function make_all(; with_examples::Bool = true) "Index" => "package_index.md", "List of Finite Elements" => "fems.md", "Base Structures" => Any[ - "fespace.md", - "fevector.md", - "fematrix.md", - "functionoperators.md", - "feevaluator.md", - "interpolations.md", - "quadrature.md" - ], + "fespace.md", + "fevector.md", + "fematrix.md", + "functionoperators.md", + "feevaluator.md", + "interpolations.md", + "quadrature.md", + ], "Advanced Stuff" => Any[ "pointevaluators.md", "segmentintegrators.md", - "plots.md" + "plots.md", ], "Tutorial Notebooks" => pluto_examples, "Examples" => module_examples, ] ) - cleanexamples() - + return cleanexamples() + end make_all(; with_examples = true) diff --git a/examples/Example200_LowLevelPoisson.jl b/examples/Example200_LowLevelPoisson.jl index d4aed552..18d08db4 100644 --- a/examples/Example200_LowLevelPoisson.jl +++ b/examples/Example200_LowLevelPoisson.jl @@ -36,176 +36,176 @@ const f = x -> x[1] - x[2] function main(; maxnref = 8, order = 2, Plotter = nothing) - ## Finite element type - FEType = H1Pk{1, 2, order} - - ## run once on a tiny mesh for compiling - X = LinRange(0, 1, 4) - xgrid = simplexgrid(X, X) - FES = FESpace{FEType}(xgrid) - sol, time_assembly, time_solve = solve_poisson_lowlevel(FES, μ, f) - - ## loop over uniform refinements + timings - plt = GridVisualizer(; Plotter = Plotter, layout = (1, 1), clear = true, resolution = (500, 500)) - loop_allocations = 0 - for level ∈ 1:maxnref - X = LinRange(0, 1, 2^level + 1) - time_grid = @elapsed xgrid = simplexgrid(X, X) - time_facenodes = @elapsed xgrid[FaceNodes] - FES = FESpace{FEType}(xgrid) - println("\nLEVEL = $level, ndofs = $(FES.ndofs)\n") - if level < 4 - println(stdout, unicode_gridplot(xgrid)) - end - time_dofmap = @elapsed FES[CellDofs] - sol, time_assembly, time_solve = solve_poisson_lowlevel(FES, μ, f) - - ## plot statistics - println(stdout, barplot(["Grid", "FaceNodes", "celldofs", "Assembly", "Solve"], [time_grid, time_facenodes, time_dofmap, time_assembly, time_solve], title = "Runtimes")) - - ## plot - if Plotter !== nothing - scalarplot!(plt[1,1], xgrid, view(sol.entries, 1:num_nodes(xgrid)), limits = (-0.0125,0.0125)) - else - sol_grad = continuify(sol[1], Gradient) - println(stdout, unicode_scalarplot(sol[1])) - end - end - - return sol, plt + ## Finite element type + FEType = H1Pk{1, 2, order} + + ## run once on a tiny mesh for compiling + X = LinRange(0, 1, 4) + xgrid = simplexgrid(X, X) + FES = FESpace{FEType}(xgrid) + sol, time_assembly, time_solve = solve_poisson_lowlevel(FES, μ, f) + + ## loop over uniform refinements + timings + plt = GridVisualizer(; Plotter = Plotter, layout = (1, 1), clear = true, resolution = (500, 500)) + loop_allocations = 0 + for level in 1:maxnref + X = LinRange(0, 1, 2^level + 1) + time_grid = @elapsed xgrid = simplexgrid(X, X) + time_facenodes = @elapsed xgrid[FaceNodes] + FES = FESpace{FEType}(xgrid) + println("\nLEVEL = $level, ndofs = $(FES.ndofs)\n") + if level < 4 + println(stdout, unicode_gridplot(xgrid)) + end + time_dofmap = @elapsed FES[CellDofs] + sol, time_assembly, time_solve = solve_poisson_lowlevel(FES, μ, f) + + ## plot statistics + println(stdout, barplot(["Grid", "FaceNodes", "celldofs", "Assembly", "Solve"], [time_grid, time_facenodes, time_dofmap, time_assembly, time_solve], title = "Runtimes")) + + ## plot + if Plotter !== nothing + scalarplot!(plt[1, 1], xgrid, view(sol.entries, 1:num_nodes(xgrid)), limits = (-0.0125, 0.0125)) + else + sol_grad = continuify(sol[1], Gradient) + println(stdout, unicode_scalarplot(sol[1])) + end + end + + return sol, plt end function solve_poisson_lowlevel(FES, μ, f) - Solution = FEVector(FES) - FES = Solution[1].FES - A = FEMatrix(FES, FES) - b = FEVector(FES) - println("Assembling...") - time_assembly = @elapsed @time begin - loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) - - ## fix boundary dofs - begin - bdofs = boundarydofs(FES) - for dof in bdofs - A.entries[dof, dof] = 1e60 - b.entries[dof] = 0 - end - end - ExtendableSparse.flush!(A.entries) - end - - ## solve - println("Solving linear system...") - time_solve = @elapsed @time copyto!(Solution.entries, A.entries \ b.entries) - - return Solution, time_assembly, time_solve + Solution = FEVector(FES) + FES = Solution[1].FES + A = FEMatrix(FES, FES) + b = FEVector(FES) + println("Assembling...") + time_assembly = @elapsed @time begin + loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) + + ## fix boundary dofs + begin + bdofs = boundarydofs(FES) + for dof in bdofs + A.entries[dof, dof] = 1.0e60 + b.entries[dof] = 0 + end + end + ExtendableSparse.flush!(A.entries) + end + + ## solve + println("Solving linear system...") + time_solve = @elapsed @time copyto!(Solution.entries, A.entries \ b.entries) + + return Solution, time_assembly, time_solve end function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) - xgrid = FES.xgrid - EG = xgrid[UniqueCellGeometries][1] - FEType = eltype(FES) - L2G = L2GTransformer(EG, xgrid, ON_CELLS) - - ## quadrature formula - qf = QuadratureRule{Float64, EG}(2 * (get_polynomialorder(FEType, EG) - 1)) - weights::Vector{Float64} = qf.w - xref::Vector{Vector{Float64}} = qf.xref - nweights::Int = length(weights) - cellvolumes = xgrid[CellVolumes] - - ## FE basis evaluator and dofmap - FEBasis_∇ = FEEvaluator(FES, Gradient, qf) - ∇vals = FEBasis_∇.cvals - FEBasis_id = FEEvaluator(FES, Identity, qf) - idvals = FEBasis_id.cvals - celldofs = FES[CellDofs] - - ## ASSEMBLY LOOP - loop_allocations = 0 - function barrier(EG, L2G::L2GTransformer) - ## barrier function to avoid allocations by type dispatch - ndofs4cell::Int = get_ndofs(ON_CELLS, FEType, EG) - Aloc = zeros(Float64, ndofs4cell, ndofs4cell) - ncells::Int = num_cells(xgrid) - dof_j::Int, dof_k::Int = 0, 0 - x::Vector{Float64} = zeros(Float64, 2) - - loop_allocations += @allocated for cell ∈ 1:ncells - ## update FE basis evaluators - FEBasis_∇.citem[] = cell - update_basis!(FEBasis_∇) - - ## assemble local stiffness matrix - for j ∈ 1:ndofs4cell, k ∈ j:ndofs4cell - temp = 0 - for qp ∈ 1:nweights - temp += weights[qp] * dot(view(∇vals, :, j, qp), view(∇vals, :, k, qp)) - end - Aloc[j, k] = temp - end - Aloc .*= μ * cellvolumes[cell] - - ## add local matrix to global matrix - for j ∈ 1:ndofs4cell - dof_j = celldofs[j, cell] - for k ∈ j:ndofs4cell - dof_k = celldofs[k, cell] - if abs(Aloc[j, k]) > 1e-15 - ## write into sparse matrix, only lines with allocations - rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k) - if k > j - rawupdateindex!(A, +, Aloc[j, k], dof_k, dof_j) - end - end - end - end - fill!(Aloc, 0) - - ## assemble right-hand side - update_trafo!(L2G, cell) - for j ∈ 1:ndofs4cell - ## right-hand side - temp = 0 - for qp ∈ 1:nweights - ## get global x for quadrature point - eval_trafo!(x, L2G, xref[qp]) - ## evaluate (f(x), v_j(x)) - temp += weights[qp] * idvals[1, j, qp] * f(x) - end - ## write into global vector - dof_j = celldofs[j, cell] - b[dof_j] += temp * cellvolumes[cell] - end - end - end - barrier(EG, L2G) - flush!(A) - return loop_allocations + xgrid = FES.xgrid + EG = xgrid[UniqueCellGeometries][1] + FEType = eltype(FES) + L2G = L2GTransformer(EG, xgrid, ON_CELLS) + + ## quadrature formula + qf = QuadratureRule{Float64, EG}(2 * (get_polynomialorder(FEType, EG) - 1)) + weights::Vector{Float64} = qf.w + xref::Vector{Vector{Float64}} = qf.xref + nweights::Int = length(weights) + cellvolumes = xgrid[CellVolumes] + + ## FE basis evaluator and dofmap + FEBasis_∇ = FEEvaluator(FES, Gradient, qf) + ∇vals = FEBasis_∇.cvals + FEBasis_id = FEEvaluator(FES, Identity, qf) + idvals = FEBasis_id.cvals + celldofs = FES[CellDofs] + + ## ASSEMBLY LOOP + loop_allocations = 0 + function barrier(EG, L2G::L2GTransformer) + ## barrier function to avoid allocations by type dispatch + ndofs4cell::Int = get_ndofs(ON_CELLS, FEType, EG) + Aloc = zeros(Float64, ndofs4cell, ndofs4cell) + ncells::Int = num_cells(xgrid) + dof_j::Int, dof_k::Int = 0, 0 + x::Vector{Float64} = zeros(Float64, 2) + + return loop_allocations += @allocated for cell in 1:ncells + ## update FE basis evaluators + FEBasis_∇.citem[] = cell + update_basis!(FEBasis_∇) + + ## assemble local stiffness matrix + for j in 1:ndofs4cell, k in j:ndofs4cell + temp = 0 + for qp in 1:nweights + temp += weights[qp] * dot(view(∇vals, :, j, qp), view(∇vals, :, k, qp)) + end + Aloc[j, k] = temp + end + Aloc .*= μ * cellvolumes[cell] + + ## add local matrix to global matrix + for j in 1:ndofs4cell + dof_j = celldofs[j, cell] + for k in j:ndofs4cell + dof_k = celldofs[k, cell] + if abs(Aloc[j, k]) > 1.0e-15 + ## write into sparse matrix, only lines with allocations + rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k) + if k > j + rawupdateindex!(A, +, Aloc[j, k], dof_k, dof_j) + end + end + end + end + fill!(Aloc, 0) + + ## assemble right-hand side + update_trafo!(L2G, cell) + for j in 1:ndofs4cell + ## right-hand side + temp = 0 + for qp in 1:nweights + ## get global x for quadrature point + eval_trafo!(x, L2G, xref[qp]) + ## evaluate (f(x), v_j(x)) + temp += weights[qp] * idvals[1, j, qp] * f(x) + end + ## write into global vector + dof_j = celldofs[j, cell] + b[dof_j] += temp * cellvolumes[cell] + end + end + end + barrier(EG, L2G) + flush!(A) + return loop_allocations end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) - ~, plt = main(; Plotter = Plotter, kwargs...) - scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example200.png"), scene; Plotter = Plotter) + ~, plt = main(; Plotter = Plotter, kwargs...) + scene = GridVisualize.reveal(plt) + return GridVisualize.save(joinpath(dir, "example200.png"), scene; Plotter = Plotter) end function runtests(; order = 2, kwargs...) #hide - FEType = H1Pk{1, 2, order} - X = LinRange(0, 1, 64) - xgrid = simplexgrid(X, X) - FES = FESpace{FEType}(xgrid) - A = FEMatrix(FES, FES) - b = FEVector(FES) - @info "ndofs = $(FES.ndofs)" - ## first assembly causes allocations when filling sparse matrix - loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) - @info "allocations in 1st assembly: $loop_allocations" - ## second assembly in same matrix should have allocation-free inner loop - loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) - @info "allocations in 2nd assembly: $loop_allocations" - @test loop_allocations == 0 + FEType = H1Pk{1, 2, order} + X = LinRange(0, 1, 64) + xgrid = simplexgrid(X, X) + FES = FESpace{FEType}(xgrid) + A = FEMatrix(FES, FES) + b = FEVector(FES) + @info "ndofs = $(FES.ndofs)" + ## first assembly causes allocations when filling sparse matrix + loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) + @info "allocations in 1st assembly: $loop_allocations" + ## second assembly in same matrix should have allocation-free inner loop + loop_allocations = assemble!(A.entries, b.entries, FES, f, μ) + @info "allocations in 2nd assembly: $loop_allocations" + return @test loop_allocations == 0 end #hide end #module diff --git a/examples/Example205_LowLevelSpaceTimePoisson.jl b/examples/Example205_LowLevelSpaceTimePoisson.jl index 1fdd59eb..5cc08ab8 100644 --- a/examples/Example205_LowLevelSpaceTimePoisson.jl +++ b/examples/Example205_LowLevelSpaceTimePoisson.jl @@ -31,57 +31,57 @@ using UnicodePlots using Test # ## data for Poisson problem -const μ = (t) -> 1e-1*t + 1*max(0,(1-2*t)) -const f = (x, t) -> sin(3*pi*x[1])*4*t - cos(3*pi*x[2])*4*(1-t) +const μ = (t) -> 1.0e-1 * t + 1 * max(0, (1 - 2 * t)) +const f = (x, t) -> sin(3 * pi * x[1]) * 4 * t - cos(3 * pi * x[2]) * 4 * (1 - t) -function main(; dt = 0.01, Tfinal = 1, level = 5, order = 1, Plotter = nothing, produce_movie = false ) +function main(; dt = 0.01, Tfinal = 1, level = 5, order = 1, Plotter = nothing, produce_movie = false) - ## Finite element type - FEType_time = H1Pk{1, 1, order} - FEType_space = H1Pk{1, 2, order} + ## Finite element type + FEType_time = H1Pk{1, 1, order} + FEType_space = H1Pk{1, 2, order} ## time grid - T = LinRange(0, Tfinal, round(Int, Tfinal/dt + 1)) + T = LinRange(0, Tfinal, round(Int, Tfinal / dt + 1)) grid_time = simplexgrid(T) ## space grid - X = LinRange(0, 1, 2^level + 1) - grid_space = simplexgrid(X, X) - + X = LinRange(0, 1, 2^level + 1) + grid_space = simplexgrid(X, X) + ## FESpaces for time and space FES_time = FESpace{FEType_time}(grid_time) FES_space = FESpace{FEType_space}(grid_space) - - ## solve + + ## solve sol = solve_poisson_lowlevel(FES_time, FES_space, μ, f) - ## visualize - if produce_movie - @info "Producing movie..." - vis=GridVisualizer(Plotter=Plotter) - movie(vis, file="example205_video.mp4") do vis - for tj = 2 : length(T) - t = T[tj] - first = (tj-1)*FES_space.ndofs+1 - last = tj*FES_space.ndofs - scalarplot!(vis, grid_space, view(sol,first:last), title = "t = $(Float16(t))") - reveal(vis) - end - end - return sol, vis - else - @info "Plotting at five times..." - plot_timesteps = [2,round(Int,length(T)/4+0.25),round(Int,length(T)/2+0.5),round(Int,length(T)-length(T)/4),FES_time.ndofs] - plt = GridVisualizer(; Plotter = Plotter, layout = (1, length(plot_timesteps)), clear = true, resolution = (200*length(plot_timesteps), 200)) - for tj = 1 : length(plot_timesteps) - t = plot_timesteps[tj] - first = (t-1)*FES_space.ndofs+1 - last = t*FES_space.ndofs - scalarplot!(plt[1,tj], grid_space, view(sol,first:last), title = "t = $(T[t])") - end - return sol, plt - end - + ## visualize + if produce_movie + @info "Producing movie..." + vis = GridVisualizer(Plotter = Plotter) + movie(vis, file = "example205_video.mp4") do vis + for tj in 2:length(T) + t = T[tj] + first = (tj - 1) * FES_space.ndofs + 1 + last = tj * FES_space.ndofs + scalarplot!(vis, grid_space, view(sol, first:last), title = "t = $(Float16(t))") + reveal(vis) + end + end + return sol, vis + else + @info "Plotting at five times..." + plot_timesteps = [2, round(Int, length(T) / 4 + 0.25), round(Int, length(T) / 2 + 0.5), round(Int, length(T) - length(T) / 4), FES_time.ndofs] + plt = GridVisualizer(; Plotter = Plotter, layout = (1, length(plot_timesteps)), clear = true, resolution = (200 * length(plot_timesteps), 200)) + for tj in 1:length(plot_timesteps) + t = plot_timesteps[tj] + first = (t - 1) * FES_space.ndofs + 1 + last = t * FES_space.ndofs + scalarplot!(plt[1, tj], grid_space, view(sol, first:last), title = "t = $(T[t])") + end + return sol, plt + end + end @@ -89,225 +89,225 @@ function solve_poisson_lowlevel(FES_time, FES_space, μ, f) ndofs_time = FES_time.ndofs ndofs_space = FES_space.ndofs ndofs_total = ndofs_time * ndofs_space - sol = zeros(Float64, ndofs_total) - - A = ExtendableSparseMatrix{Float64, Int64}(ndofs_total, ndofs_total) - b = zeros(Float64, ndofs_total) - - println("Assembling...") - time_assembly = @elapsed @time begin - loop_allocations = assemble!(A, b, FES_time, FES_space, f, μ) - - ## fix homogeneous boundary dofs - bdofs = boundarydofs(FES_space) - for sdof in bdofs - for dof_t = 1 : ndofs_time - dof = (dof_t-1)*ndofs_space + sdof - A[dof, dof] = 1e60 - b[dof] = 0 - end - end + sol = zeros(Float64, ndofs_total) + + A = ExtendableSparseMatrix{Float64, Int64}(ndofs_total, ndofs_total) + b = zeros(Float64, ndofs_total) + + println("Assembling...") + time_assembly = @elapsed @time begin + loop_allocations = assemble!(A, b, FES_time, FES_space, f, μ) + + ## fix homogeneous boundary dofs + bdofs = boundarydofs(FES_space) + for sdof in bdofs + for dof_t in 1:ndofs_time + dof = (dof_t - 1) * ndofs_space + sdof + A[dof, dof] = 1.0e60 + b[dof] = 0 + end + end ## fix initial value by zero - for j=1:ndofs_space - A[j,j] = 1e60 - b[j] = 0 - end - ExtendableSparse.flush!(A) - end + for j in 1:ndofs_space + A[j, j] = 1.0e60 + b[j] = 0 + end + ExtendableSparse.flush!(A) + end @info ".... spy plot of system matrix:\n$(UnicodePlots.spy(sparse(A.cscmatrix)))" - ## solve - println("Solving linear system...") - time_solve = @elapsed @time copyto!(sol, A \ b) + ## solve + println("Solving linear system...") + time_solve = @elapsed @time copyto!(sol, A \ b) - ## compute linear residual - @show norm(A*sol - b) + ## compute linear residual + @show norm(A * sol - b) - return sol + return sol end function assemble!(A::ExtendableSparseMatrix, b::Vector, FES_time, FES_space, f, μ = 1) - - ## get space and time grids - grid_time = FES_time.xgrid - grid_space = FES_space.xgrid - ## get number of degrees of freedom + ## get space and time grids + grid_time = FES_time.xgrid + grid_space = FES_space.xgrid + + ## get number of degrees of freedom ndofs_time = FES_time.ndofs ndofs_space = FES_space.ndofs - ## get local to global maps - EG_time = grid_time[UniqueCellGeometries][1] - EG_space = grid_space[UniqueCellGeometries][1] - L2G_time = L2GTransformer(EG_time, grid_time, ON_CELLS) - L2G_space = L2GTransformer(EG_space, grid_space, ON_CELLS) - - ## get finite element types - FEType_time = eltype(FES_time) - FEType_space = eltype(FES_space) - - ## quadrature formula in space - qf_space = QuadratureRule{Float64, EG_space}(2 * (get_polynomialorder(FEType_space, EG_space) - 1)) - weights_space::Vector{Float64} = qf_space.w - xref_space::Vector{Vector{Float64}} = qf_space.xref - nweights_space::Int = length(weights_space) - cellvolumes_space = grid_space[CellVolumes] - - ## quadrature formula in time - qf_time = QuadratureRule{Float64, EG_time}(2 * (get_polynomialorder(FEType_time, EG_time) - 1)) - weights_time::Vector{Float64} = qf_time.w - xref_time::Vector{Vector{Float64}} = qf_time.xref - nweights_time::Int = length(weights_time) - cellvolumes_time = grid_time[CellVolumes] - - ## FE basis evaluators and dofmap for space elements - FEBasis_space_∇ = FEEvaluator(FES_space, Gradient, qf_space) - ∇vals_space = FEBasis_space_∇.cvals - FEBasis_space_id = FEEvaluator(FES_space, Identity, qf_space) - idvals_space = FEBasis_space_id.cvals - celldofs_space = FES_space[ExtendableFEMBase.CellDofs] - - ## FE basis evaluators and dofmap for time elements - FEBasis_time_∇ = FEEvaluator(FES_time, Gradient, qf_time) - ∇vals_time = FEBasis_time_∇.cvals - FEBasis_time_id = FEEvaluator(FES_time, Identity, qf_time) - idvals_time = FEBasis_time_id.cvals - celldofs_time = FES_time[ExtendableFEMBase.CellDofs] - - ## ASSEMBLY LOOP - loop_allocations = 0 - function barrier(EG_time, EG_space, L2G_time::L2GTransformer, L2G_space::L2GTransformer) - ## barrier function to avoid allocations by type dispatch - - ndofs4cell_time::Int = get_ndofs(ON_CELLS, FEType_time, EG_time) - ndofs4cell_space::Int = get_ndofs(ON_CELLS, FEType_space, EG_space) - Aloc = zeros(Float64, ndofs4cell_space, ndofs4cell_space) + ## get local to global maps + EG_time = grid_time[UniqueCellGeometries][1] + EG_space = grid_space[UniqueCellGeometries][1] + L2G_time = L2GTransformer(EG_time, grid_time, ON_CELLS) + L2G_space = L2GTransformer(EG_space, grid_space, ON_CELLS) + + ## get finite element types + FEType_time = eltype(FES_time) + FEType_space = eltype(FES_space) + + ## quadrature formula in space + qf_space = QuadratureRule{Float64, EG_space}(2 * (get_polynomialorder(FEType_space, EG_space) - 1)) + weights_space::Vector{Float64} = qf_space.w + xref_space::Vector{Vector{Float64}} = qf_space.xref + nweights_space::Int = length(weights_space) + cellvolumes_space = grid_space[CellVolumes] + + ## quadrature formula in time + qf_time = QuadratureRule{Float64, EG_time}(2 * (get_polynomialorder(FEType_time, EG_time) - 1)) + weights_time::Vector{Float64} = qf_time.w + xref_time::Vector{Vector{Float64}} = qf_time.xref + nweights_time::Int = length(weights_time) + cellvolumes_time = grid_time[CellVolumes] + + ## FE basis evaluators and dofmap for space elements + FEBasis_space_∇ = FEEvaluator(FES_space, Gradient, qf_space) + ∇vals_space = FEBasis_space_∇.cvals + FEBasis_space_id = FEEvaluator(FES_space, Identity, qf_space) + idvals_space = FEBasis_space_id.cvals + celldofs_space = FES_space[ExtendableFEMBase.CellDofs] + + ## FE basis evaluators and dofmap for time elements + FEBasis_time_∇ = FEEvaluator(FES_time, Gradient, qf_time) + ∇vals_time = FEBasis_time_∇.cvals + FEBasis_time_id = FEEvaluator(FES_time, Identity, qf_time) + idvals_time = FEBasis_time_id.cvals + celldofs_time = FES_time[ExtendableFEMBase.CellDofs] + + ## ASSEMBLY LOOP + loop_allocations = 0 + function barrier(EG_time, EG_space, L2G_time::L2GTransformer, L2G_space::L2GTransformer) + ## barrier function to avoid allocations by type dispatch + + ndofs4cell_time::Int = get_ndofs(ON_CELLS, FEType_time, EG_time) + ndofs4cell_space::Int = get_ndofs(ON_CELLS, FEType_space, EG_space) + Aloc = zeros(Float64, ndofs4cell_space, ndofs4cell_space) Mloc = zeros(Float64, ndofs4cell_time, ndofs4cell_time) - ncells_space::Int = num_cells(grid_space) + ncells_space::Int = num_cells(grid_space) ncells_time::Int = num_cells(grid_time) - x::Vector{Float64} = zeros(Float64, 2) - t::Vector{Float64} = zeros(Float64, 1) + x::Vector{Float64} = zeros(Float64, 2) + t::Vector{Float64} = zeros(Float64, 1) ## assemble Laplacian - loop_allocations += @allocated for cell ∈ 1:ncells_space - ## update FE basis evaluators for space - FEBasis_space_∇.citem[] = cell - update_basis!(FEBasis_space_∇) - - ## assemble local stiffness matrix in space - for j ∈ 1:ndofs4cell_space, k ∈ 1:ndofs4cell_space - temp = 0 - for qp ∈ 1:nweights_space - temp += weights_space[qp] * dot(view(∇vals_space, :, j, qp), view(∇vals_space, :, k, qp)) - end - Aloc[j, k] = temp - end - Aloc .*= cellvolumes_space[cell] - - ## add local matrix to global matrix - for time_cell ∈ 1:ncells_time - update_trafo!(L2G_time, time_cell) - for jT ∈ 1:ndofs4cell_time, kT ∈ 1:ndofs4cell_time + loop_allocations += @allocated for cell in 1:ncells_space + ## update FE basis evaluators for space + FEBasis_space_∇.citem[] = cell + update_basis!(FEBasis_space_∇) + + ## assemble local stiffness matrix in space + for j in 1:ndofs4cell_space, k in 1:ndofs4cell_space + temp = 0 + for qp in 1:nweights_space + temp += weights_space[qp] * dot(view(∇vals_space, :, j, qp), view(∇vals_space, :, k, qp)) + end + Aloc[j, k] = temp + end + Aloc .*= cellvolumes_space[cell] + + ## add local matrix to global matrix + for time_cell in 1:ncells_time + update_trafo!(L2G_time, time_cell) + for jT in 1:ndofs4cell_time, kT in 1:ndofs4cell_time dofTj = celldofs_time[jT, time_cell] dofTk = celldofs_time[kT, time_cell] - for qpT ∈ 1:nweights_time - ## evaluate time coordinate and μ - eval_trafo!(t, L2G_time, xref_time[qpT]) - factor = μ(t[1]) * weights_time[qpT] * idvals_time[1,jT,qpT] * idvals_time[1,kT,qpT] * cellvolumes_time[time_cell] - for j ∈ 1:ndofs4cell_space + for qpT in 1:nweights_time + ## evaluate time coordinate and μ + eval_trafo!(t, L2G_time, xref_time[qpT]) + factor = μ(t[1]) * weights_time[qpT] * idvals_time[1, jT, qpT] * idvals_time[1, kT, qpT] * cellvolumes_time[time_cell] + for j in 1:ndofs4cell_space dof_j = celldofs_space[j, cell] + (dofTj - 1) * ndofs_space - for k ∈ 1:ndofs4cell_space + for k in 1:ndofs4cell_space dof_k = celldofs_space[k, cell] + (dofTk - 1) * ndofs_space - if abs(Aloc[j, k]) > 1e-15 + if abs(Aloc[j, k]) > 1.0e-15 ## write into sparse matrix, only lines with allocations - rawupdateindex!(A, +, Aloc[j, k]*factor, dof_j, dof_k) + rawupdateindex!(A, +, Aloc[j, k] * factor, dof_j, dof_k) end end end - end + end end end - fill!(Aloc, 0) - - ## assemble right-hand side - update_trafo!(L2G_space, cell) - for qp ∈ 1:nweights_space - ## evaluate coordinates of quadrature point in space - eval_trafo!(x, L2G_space, xref_space[qp]) - for time_cell ∈ 1:ncells_time - update_trafo!(L2G_time, time_cell) - for qpT ∈ 1:nweights_time - ## evaluate time coordinate - eval_trafo!(t, L2G_time, xref_time[qpT]) - - ## evaluate right-hand side in x and t - fval = f(x, t[1]) - - ## multiply with test function and add to right-hand side - for j ∈ 1:ndofs4cell_space - temp = weights_time[qpT] * weights_space[qp] * idvals_space[1, j, qp] * fval * cellvolumes_space[cell] * cellvolumes_time[time_cell] - - ## write into global vector - for jT ∈ 1:ndofs4cell_time - dof_j = celldofs_space[j, cell] + (celldofs_time[jT, time_cell] - 1) * ndofs_space - b[dof_j] += temp * idvals_time[1, jT, qpT] - end - end - end + fill!(Aloc, 0) + + ## assemble right-hand side + 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 in 1:ncells_time + update_trafo!(L2G_time, time_cell) + for qpT in 1:nweights_time + ## evaluate time coordinate + eval_trafo!(t, L2G_time, xref_time[qpT]) + + ## evaluate right-hand side in x and t + fval = f(x, t[1]) + + ## multiply with test function and add to right-hand side + for j in 1:ndofs4cell_space + temp = weights_time[qpT] * weights_space[qp] * idvals_space[1, j, qp] * fval * cellvolumes_space[cell] * cellvolumes_time[time_cell] + + ## write into global vector + for jT in 1:ndofs4cell_time + dof_j = celldofs_space[j, cell] + (celldofs_time[jT, time_cell] - 1) * ndofs_space + b[dof_j] += temp * idvals_time[1, jT, qpT] + end + end + end end - end - end + end + end ## assemble time derivative - loop_allocations += @allocated for time_cell ∈ 1:ncells_time - ## update FE basis evaluators for time derivative - FEBasis_time_∇.citem[] = time_cell - update_basis!(FEBasis_time_∇) - - ## assemble local convection term in time - for j ∈ 1:ndofs4cell_time, k ∈ 1:ndofs4cell_time - temp = 0 - for qpT ∈ 1:nweights_time - temp += weights_time[qpT] * dot(view(∇vals_time, :, j, qpT), view(∇vals_time, :, k, qpT)) - end - Mloc[j, k] = temp - end - Mloc .*= cellvolumes_time[time_cell] - - ## add local matrix to global matrix - for cell ∈ 1:ncells_space - for jX ∈ 1:ndofs4cell_space, kX ∈ 1:ndofs4cell_space + return loop_allocations += @allocated for time_cell in 1:ncells_time + ## update FE basis evaluators for time derivative + FEBasis_time_∇.citem[] = time_cell + update_basis!(FEBasis_time_∇) + + ## assemble local convection term in time + for j in 1:ndofs4cell_time, k in 1:ndofs4cell_time + temp = 0 + for qpT in 1:nweights_time + temp += weights_time[qpT] * dot(view(∇vals_time, :, j, qpT), view(∇vals_time, :, k, qpT)) + end + Mloc[j, k] = temp + end + Mloc .*= cellvolumes_time[time_cell] + + ## add local matrix to global matrix + for cell in 1:ncells_space + for jX in 1:ndofs4cell_space, kX in 1:ndofs4cell_space dofXj = celldofs_space[jX, cell] dofXk = celldofs_space[kX, cell] - for qpX ∈ 1:nweights_space + for qpX in 1:nweights_space factor = weights_space[qpX] * idvals_space[1, jX, qpX] * idvals_space[1, kX, qpX] * cellvolumes_space[cell] - for j ∈ 1:ndofs4cell_time + for j in 1:ndofs4cell_time dof_j = dofXj + (celldofs_time[j, time_cell] - 1) * ndofs_space - for k ∈ 1:ndofs4cell_time + for k in 1:ndofs4cell_time dof_k = dofXk + (celldofs_time[k, time_cell] - 1) * ndofs_space - if abs(Mloc[j, k]) > 1e-15 + if abs(Mloc[j, k]) > 1.0e-15 ## write into sparse matrix, only lines with allocations - rawupdateindex!(A, +, Mloc[j, k]*factor, dof_j, dof_k) + rawupdateindex!(A, +, Mloc[j, k] * factor, dof_j, dof_k) end end end - end + end end end - fill!(Mloc, 0) - end - end - barrier(EG_time, EG_space, L2G_time, L2G_space) - flush!(A) - return loop_allocations + fill!(Mloc, 0) + end + end + barrier(EG_time, EG_space, L2G_time, L2G_space) + flush!(A) + return loop_allocations end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) - ~, plt = main(; Plotter = Plotter, kwargs...) - scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example205.png"), scene; Plotter = Plotter) + ~, plt = main(; Plotter = Plotter, kwargs...) + scene = GridVisualize.reveal(plt) + return GridVisualize.save(joinpath(dir, "example205.png"), scene; Plotter = Plotter) end end #module diff --git a/examples/Example210_LowLevelNavierStokes.jl b/examples/Example210_LowLevelNavierStokes.jl index 7153a298..48a37f6e 100644 --- a/examples/Example210_LowLevelNavierStokes.jl +++ b/examples/Example210_LowLevelNavierStokes.jl @@ -39,356 +39,360 @@ using ForwardDiff using DiffResults ## data for Poisson problem -const μ = 1e-2 +const μ = 1.0e-2 function f!(fval, x, t) # right-hand side - fval[1] = 8.0 * π * π * μ * exp(-8.0 * π * π * μ * t) * sin(2.0 * π * x[1]) * sin(2.0 * π * x[2]) - fval[2] = 8.0 * π * π * μ * exp(-8.0 * π * π * μ * t) * cos(2.0 * π * x[1]) * cos(2.0 * π * x[2]) - return nothing + fval[1] = 8.0 * π * π * μ * exp(-8.0 * π * π * μ * t) * sin(2.0 * π * x[1]) * sin(2.0 * π * x[2]) + fval[2] = 8.0 * π * π * μ * exp(-8.0 * π * π * μ * t) * cos(2.0 * π * x[1]) * cos(2.0 * π * x[2]) + return nothing end ## exact velocity (for boundary data and error calculation) function u!(uval, qpinfo) - x = qpinfo.x - t = qpinfo.time - uval[1] = exp(-8.0 * π * π * μ * t) * sin(2.0 * π * x[1]) * sin(2.0 * π * x[2]) - uval[2] = exp(-8.0 * π * π * μ * t) * cos(2.0 * π * x[1]) * cos(2.0 * π * x[2]) - return nothing + x = qpinfo.x + t = qpinfo.time + uval[1] = exp(-8.0 * π * π * μ * t) * sin(2.0 * π * x[1]) * sin(2.0 * π * x[2]) + uval[2] = exp(-8.0 * π * π * μ * t) * cos(2.0 * π * x[1]) * cos(2.0 * π * x[2]) + return nothing end ## exact pressure (for error calculation) function p!(pval, qpinfo) - x = qpinfo.x - t = qpinfo.time - pval[1] = exp(-16 * pi * pi * μ * t) * (cos(4 * pi * x[1]) - cos(4 * pi * x[2])) / 4 - return nothing + x = qpinfo.x + t = qpinfo.time + pval[1] = exp(-16 * pi * pi * μ * t) * (cos(4 * pi * x[1]) - cos(4 * pi * x[2])) / 4 + return nothing end function main(; nref = 5, teval = 0, order = 2, Plotter = nothing) - @assert order >= 2 - - ## create grid - X = LinRange(0, 1, 2^nref + 1) - Y = LinRange(0, 1, 2^nref + 1) - println("Creating grid...") - @time xgrid = simplexgrid(X, Y) - - ## create FESpace - println("Creating FESpace...") - FETypes = [H1Pk{2, 2, order}, H1Pk{1, 2, order - 1}] - @time FES = [FESpace{FETypes[1]}(xgrid; name = "velocity space"), - FESpace{FETypes[2]}(xgrid; name = "pressure space")] - FES - - ## solve - sol = solve_stokes_lowlevel(FES; teval = teval) - - ## move integral mean of pressure - pmean = sum(compute_error(sol[2], nothing, order, 1)) - for j ∈ 1:sol[2].FES.ndofs - sol[2][j] -= pmean - end - - ## 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)") - println("l2 error pressure = $(error_p)") - - ## plot - plt = GridVisualizer(; Plotter = Plotter, layout = (1, 1), clear = true, resolution = (500, 500)) - scalarplot!(plt[1, 1], xgrid, nodevalues(sol[1]; abs = true)[1, :]; title = "|u| + quiver", Plotter = Plotter) - vectorplot!(plt[1, 1], xgrid, eval_func_bary(PointEvaluator([(1, Identity)], sol)), clear = false) - - return sol, plt + @assert order >= 2 + + ## create grid + X = LinRange(0, 1, 2^nref + 1) + Y = LinRange(0, 1, 2^nref + 1) + println("Creating grid...") + @time xgrid = simplexgrid(X, Y) + + ## create FESpace + println("Creating FESpace...") + FETypes = [H1Pk{2, 2, order}, H1Pk{1, 2, order - 1}] + @time FES = [ + FESpace{FETypes[1]}(xgrid; name = "velocity space"), + FESpace{FETypes[2]}(xgrid; name = "pressure space"), + ] + FES + + ## solve + 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 + + ## 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)") + println("l2 error pressure = $(error_p)") + + ## plot + plt = GridVisualizer(; Plotter = Plotter, layout = (1, 1), clear = true, resolution = (500, 500)) + scalarplot!(plt[1, 1], xgrid, nodevalues(sol[1]; abs = true)[1, :]; title = "|u| + quiver", Plotter = Plotter) + vectorplot!(plt[1, 1], xgrid, eval_func_bary(PointEvaluator([(1, Identity)], sol)), clear = false) + + return sol, plt end ## computes error and integrals function compute_error(uh::FEVectorBlock, u, order = get_polynomialorder(get_FEType(uh), uh.FES.xgrid[CellGeometries][1]), p = 2) - xgrid = uh.FES.xgrid - FES = uh.FES - EG = xgrid[UniqueCellGeometries][1] - ncomponents = get_ncomponents(uh) - cellvolumes = xgrid[CellVolumes] - celldofs = FES[CellDofs] - error = zeros(Float64, ncomponents, num_cells(xgrid)) - uhval = zeros(Float64, ncomponents) - uval = zeros(Float64, ncomponents) - L2G = L2GTransformer(EG, xgrid, ON_CELLS) - QP = QPInfos(xgrid) - qf = VertexRule(EG, order) - FEB = FEEvaluator(FES, Identity, qf) - - function barrier(L2G::L2GTransformer) - for cell in 1:num_cells(xgrid) - update_trafo!(L2G, cell) - update_basis!(FEB, cell) - for (qp, weight) in enumerate(qf.w) - ## evaluate uh - fill!(uhval, 0) - eval_febe!(uhval, FEB, view(view(uh), view(celldofs, :, cell)), qp) - - ## evaluate u - if u !== nothing - fill!(uval, 0) - eval_trafo!(QP.x, L2G, qf.xref[qp]) - u(uval, QP) - end - - ## evaluate error - for d ∈ 1:ncomponents - error[d, cell] += (uhval[d] - uval[d]) .^ p * cellvolumes[cell] * weight - end - end - end - end - - barrier(L2G) - return error + xgrid = uh.FES.xgrid + FES = uh.FES + EG = xgrid[UniqueCellGeometries][1] + ncomponents = get_ncomponents(uh) + cellvolumes = xgrid[CellVolumes] + celldofs = FES[CellDofs] + error = zeros(Float64, ncomponents, num_cells(xgrid)) + uhval = zeros(Float64, ncomponents) + uval = zeros(Float64, ncomponents) + L2G = L2GTransformer(EG, xgrid, ON_CELLS) + QP = QPInfos(xgrid) + qf = VertexRule(EG, order) + FEB = FEEvaluator(FES, Identity, qf) + + function barrier(L2G::L2GTransformer) + for cell in 1:num_cells(xgrid) + update_trafo!(L2G, cell) + update_basis!(FEB, cell) + for (qp, weight) in enumerate(qf.w) + ## evaluate uh + fill!(uhval, 0) + eval_febe!(uhval, FEB, view(view(uh), view(celldofs, :, cell)), qp) + + ## evaluate u + if u !== nothing + fill!(uval, 0) + eval_trafo!(QP.x, L2G, qf.xref[qp]) + u(uval, QP) + end + + ## evaluate error + for d in 1:ncomponents + error[d, cell] += (uhval[d] - uval[d]) .^ p * cellvolumes[cell] * weight + end + end + end + return + end + + barrier(L2G) + return error end function solve_stokes_lowlevel(FES; teval = 0) - println("Initializing system...") - sol = FEVector(FES) - A = FEMatrix(FES) - b = FEVector(FES) - @time update_system! = prepare_assembly!(A, b, FES[1], FES[2], sol) - @time update_system!(true, false) - Alin = deepcopy(A) ## = keep linear part of system matrix - blin = deepcopy(b) ## = keep linear part of right-hand side - - println("Prepare boundary conditions...") - @time begin - u_init = FEVector(FES) - interpolate!(u_init[1], u!; time = teval) - fixed_dofs = boundarydofs(FES[1]) - push!(fixed_dofs, FES[1].ndofs + 1) ## fix one pressure dof - end - - for it ∈ 1:20 - ## solve - println("\nITERATION $it\n=============") - println("Solving linear system...") - @time copyto!(sol.entries, A.entries \ b.entries) - res = A.entries.cscmatrix * sol.entries .- b.entries - for dof in fixed_dofs - res[dof] = 0 - end - linres = norm(res) - println("linear residual = $linres") - - fill!(A.entries.cscmatrix.nzval, 0) - fill!(b.entries, 0) - println("Updating linear system...") - @time begin - update_system!(false, true) - A.entries.cscmatrix += Alin.entries.cscmatrix - b.entries .+= blin.entries - end - - ## fix boundary dofs - for dof in fixed_dofs - A.entries[dof, dof] = 1e60 - b.entries[dof] = 1e60 * u_init.entries[dof] - end - ExtendableSparse.flush!(A.entries) - - ## calculate nonlinear residual - res = A.entries.cscmatrix * sol.entries .- b.entries - for dof in fixed_dofs - res[dof] = 0 - end - nlres = norm(res) - println("nonlinear residual = $nlres") - if nlres < max(1e-12, 20 * linres) - break - end - end - - return sol + println("Initializing system...") + sol = FEVector(FES) + A = FEMatrix(FES) + b = FEVector(FES) + @time update_system! = prepare_assembly!(A, b, FES[1], FES[2], sol) + @time update_system!(true, false) + Alin = deepcopy(A) ## = keep linear part of system matrix + blin = deepcopy(b) ## = keep linear part of right-hand side + + println("Prepare boundary conditions...") + @time begin + u_init = FEVector(FES) + interpolate!(u_init[1], u!; time = teval) + fixed_dofs = boundarydofs(FES[1]) + push!(fixed_dofs, FES[1].ndofs + 1) ## fix one pressure dof + end + + for it in 1:20 + ## solve + println("\nITERATION $it\n=============") + println("Solving linear system...") + @time copyto!(sol.entries, A.entries \ b.entries) + res = A.entries.cscmatrix * sol.entries .- b.entries + for dof in fixed_dofs + res[dof] = 0 + end + linres = norm(res) + println("linear residual = $linres") + + fill!(A.entries.cscmatrix.nzval, 0) + fill!(b.entries, 0) + println("Updating linear system...") + @time begin + update_system!(false, true) + A.entries.cscmatrix += Alin.entries.cscmatrix + b.entries .+= blin.entries + 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) + + ## calculate nonlinear residual + res = A.entries.cscmatrix * sol.entries .- b.entries + for dof in fixed_dofs + res[dof] = 0 + end + nlres = norm(res) + println("nonlinear residual = $nlres") + if nlres < max(1.0e-12, 20 * linres) + break + end + end + + return sol end function prepare_assembly!(A, b, FESu, FESp, sol; teval = 0) - A = A.entries - b = b.entries - sol = sol.entries - xgrid = FESu.xgrid - EG = xgrid[UniqueCellGeometries][1] - FEType_u = eltype(FESu) - FEType_p = eltype(FESp) - L2G = L2GTransformer(EG, xgrid, ON_CELLS) - cellvolumes = xgrid[CellVolumes] - ncells::Int = num_cells(xgrid) - - ## dofmap - CellDofs_u = FESu[ExtendableFEMBase.CellDofs] - CellDofs_p = FESp[ExtendableFEMBase.CellDofs] - offset_p = FESu.ndofs - - ## quadrature formula - qf = QuadratureRule{Float64, EG}(3 * get_polynomialorder(FEType_u, EG) - 1) - weights::Vector{Float64} = qf.w - xref::Vector{Vector{Float64}} = qf.xref - nweights::Int = length(weights) - - ## FE basis evaluator - FEBasis_∇u = FEEvaluator(FESu, Gradient, qf) - ∇uvals = FEBasis_∇u.cvals - FEBasis_idu = FEEvaluator(FESu, Identity, qf) - iduvals = FEBasis_idu.cvals - FEBasis_idp = FEEvaluator(FESp, Identity, qf) - idpvals = FEBasis_idp.cvals - - ## prepare automatic differentation of convection operator - function operator!(result, input) - ## result = (u ⋅ ∇)u - result[1] = input[1] * input[3] + input[2] * input[4] - result[2] = input[1] * input[5] + input[2] * input[6] - end - result = Vector{Float64}(undef, 2) - input = Vector{Float64}(undef, 6) - tempV = zeros(Float64, 2) - Dresult = DiffResults.JacobianResult(result, input) - cfg = ForwardDiff.JacobianConfig(operator!, result, input, ForwardDiff.Chunk{6}()) - jac = DiffResults.jacobian(Dresult) - value = DiffResults.value(Dresult) - - ## ASSEMBLY LOOP - function barrier(EG, L2G::L2GTransformer, linear::Bool, nonlinear::Bool) - ## barrier function to avoid allocations caused by L2G - - ndofs4cell_u::Int = get_ndofs(ON_CELLS, FEType_u, EG) - ndofs4cell_p::Int = get_ndofs(ON_CELLS, FEType_p, EG) - Aloc = zeros(Float64, ndofs4cell_u, ndofs4cell_u) - Bloc = zeros(Float64, ndofs4cell_u, ndofs4cell_p) - dof_j::Int, dof_k::Int = 0, 0 - fval::Vector{Float64} = zeros(Float64, 2) - x::Vector{Float64} = zeros(Float64, 2) - - for cell ∈ 1:ncells - ## update FE basis evaluators - update_basis!(FEBasis_∇u, cell) - update_basis!(FEBasis_idu, cell) - update_basis!(FEBasis_idp, cell) - - ## assemble local stiffness matrix (symmetric) - if (linear) - for j ∈ 1:ndofs4cell_u, k ∈ 1:ndofs4cell_u - temp = 0 - for qp ∈ 1:nweights - temp += weights[qp] * dot(view(∇uvals, :, j, qp), view(∇uvals, :, k, qp)) - end - Aloc[k, j] = μ * temp - end - - ## assemble div-pressure coupling - for j ∈ 1:ndofs4cell_u, k ∈ 1:ndofs4cell_p - temp = 0 - for qp ∈ 1:nweights - temp -= weights[qp] * (∇uvals[1, j, qp] + ∇uvals[4, j, qp]) * - idpvals[1, k, qp] - end - Bloc[j, k] = temp - end - Bloc .*= cellvolumes[cell] - - ## assemble right-hand side - update_trafo!(L2G, cell) - for j ∈ 1:ndofs4cell_u - ## right-hand side - temp = 0 - for qp ∈ 1:nweights - ## get global x for quadrature point - eval_trafo!(x, L2G, xref[qp]) - ## evaluate (f(x), v_j(x)) - f!(fval, x, teval) - temp += weights[qp] * dot(view(iduvals, :, j, qp), fval) - end - ## write into global vector - dof_j = CellDofs_u[j, cell] - b[dof_j] += temp * cellvolumes[cell] - end - end - - ## assemble nonlinear term - if (nonlinear) - for qp ∈ 1:nweights - fill!(input, 0) - for j ∈ 1:ndofs4cell_u - dof_j = CellDofs_u[j, cell] - for d ∈ 1:2 - input[d] += sol[dof_j] * iduvals[d, j, qp] - end - for d ∈ 1:4 - input[2+d] += sol[dof_j] * ∇uvals[d, j, qp] - end - end - - ## evaluate jacobian - ForwardDiff.chunk_mode_jacobian!(Dresult, operator!, result, input, cfg) - - ## update matrix - for j ∈ 1:ndofs4cell_u - ## multiply ansatz function with local jacobian - fill!(tempV, 0) - for d ∈ 1:2 - tempV[1] += jac[1, d] * iduvals[d, j, qp] - tempV[2] += jac[2, d] * iduvals[d, j, qp] - end - for d ∈ 1:4 - tempV[1] += jac[1, 2+d] * ∇uvals[d, j, qp] - tempV[2] += jac[2, 2+d] * ∇uvals[d, j, qp] - end - - ## multiply test function operator evaluation - for k ∈ 1:ndofs4cell_u - Aloc[k, j] += dot(tempV, view(iduvals, :, k, qp)) * weights[qp] - end - end - - ## update rhs - mul!(tempV, jac, input) - tempV .-= value - for j ∈ 1:ndofs4cell_u - dof_j = CellDofs_u[j, cell] - b[dof_j] += dot(tempV, view(iduvals, :, j, qp)) * weights[qp] * cellvolumes[cell] - end - end - end - - ## add local matrices to global matrix - Aloc .*= cellvolumes[cell] - for j ∈ 1:ndofs4cell_u - dof_j = CellDofs_u[j, cell] - for k ∈ 1:ndofs4cell_u - dof_k = CellDofs_u[k, cell] - rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k) - end - if (linear) - for k ∈ 1:ndofs4cell_p - dof_k = CellDofs_p[k, cell] + offset_p - rawupdateindex!(A, +, Bloc[j, k], dof_j, dof_k) - rawupdateindex!(A, +, Bloc[j, k], dof_k, dof_j) - end - end - end - fill!(Aloc, 0) - fill!(Bloc, 0) - end - end - - function update_system!(linear::Bool, nonlinear::Bool) - barrier(EG, L2G, linear, nonlinear) - flush!(A) - end - update_system! + A = A.entries + b = b.entries + sol = sol.entries + xgrid = FESu.xgrid + EG = xgrid[UniqueCellGeometries][1] + FEType_u = eltype(FESu) + FEType_p = eltype(FESp) + L2G = L2GTransformer(EG, xgrid, ON_CELLS) + cellvolumes = xgrid[CellVolumes] + ncells::Int = num_cells(xgrid) + + ## dofmap + CellDofs_u = FESu[ExtendableFEMBase.CellDofs] + CellDofs_p = FESp[ExtendableFEMBase.CellDofs] + offset_p = FESu.ndofs + + ## quadrature formula + qf = QuadratureRule{Float64, EG}(3 * get_polynomialorder(FEType_u, EG) - 1) + weights::Vector{Float64} = qf.w + xref::Vector{Vector{Float64}} = qf.xref + nweights::Int = length(weights) + + ## FE basis evaluator + FEBasis_∇u = FEEvaluator(FESu, Gradient, qf) + ∇uvals = FEBasis_∇u.cvals + FEBasis_idu = FEEvaluator(FESu, Identity, qf) + iduvals = FEBasis_idu.cvals + FEBasis_idp = FEEvaluator(FESp, Identity, qf) + idpvals = FEBasis_idp.cvals + + ## prepare automatic differentation of convection operator + function operator!(result, input) + ## result = (u ⋅ ∇)u + result[1] = input[1] * input[3] + input[2] * input[4] + return result[2] = input[1] * input[5] + input[2] * input[6] + end + result = Vector{Float64}(undef, 2) + input = Vector{Float64}(undef, 6) + tempV = zeros(Float64, 2) + Dresult = DiffResults.JacobianResult(result, input) + cfg = ForwardDiff.JacobianConfig(operator!, result, input, ForwardDiff.Chunk{6}()) + jac = DiffResults.jacobian(Dresult) + value = DiffResults.value(Dresult) + + ## ASSEMBLY LOOP + function barrier(EG, L2G::L2GTransformer, linear::Bool, nonlinear::Bool) + ## barrier function to avoid allocations caused by L2G + + ndofs4cell_u::Int = get_ndofs(ON_CELLS, FEType_u, EG) + ndofs4cell_p::Int = get_ndofs(ON_CELLS, FEType_p, EG) + Aloc = zeros(Float64, ndofs4cell_u, ndofs4cell_u) + Bloc = zeros(Float64, ndofs4cell_u, ndofs4cell_p) + dof_j::Int, dof_k::Int = 0, 0 + fval::Vector{Float64} = zeros(Float64, 2) + x::Vector{Float64} = zeros(Float64, 2) + + for cell in 1:ncells + ## update FE basis evaluators + update_basis!(FEBasis_∇u, cell) + update_basis!(FEBasis_idu, cell) + update_basis!(FEBasis_idp, cell) + + ## assemble local stiffness matrix (symmetric) + if (linear) + for j in 1:ndofs4cell_u, k in 1:ndofs4cell_u + temp = 0 + for qp in 1:nweights + temp += weights[qp] * dot(view(∇uvals, :, j, qp), view(∇uvals, :, k, qp)) + end + Aloc[k, j] = μ * temp + end + + ## assemble div-pressure coupling + for j in 1:ndofs4cell_u, k in 1:ndofs4cell_p + temp = 0 + for qp in 1:nweights + temp -= weights[qp] * (∇uvals[1, j, qp] + ∇uvals[4, j, qp]) * + idpvals[1, k, qp] + end + Bloc[j, k] = temp + end + Bloc .*= cellvolumes[cell] + + ## assemble right-hand side + update_trafo!(L2G, cell) + for j in 1:ndofs4cell_u + ## right-hand side + temp = 0 + for qp in 1:nweights + ## get global x for quadrature point + eval_trafo!(x, L2G, xref[qp]) + ## evaluate (f(x), v_j(x)) + f!(fval, x, teval) + temp += weights[qp] * dot(view(iduvals, :, j, qp), fval) + end + ## write into global vector + dof_j = CellDofs_u[j, cell] + b[dof_j] += temp * cellvolumes[cell] + end + end + + ## assemble nonlinear term + if (nonlinear) + for qp in 1:nweights + fill!(input, 0) + for j in 1:ndofs4cell_u + dof_j = CellDofs_u[j, cell] + for d in 1:2 + input[d] += sol[dof_j] * iduvals[d, j, qp] + end + for d in 1:4 + input[2 + d] += sol[dof_j] * ∇uvals[d, j, qp] + end + end + + ## evaluate jacobian + ForwardDiff.chunk_mode_jacobian!(Dresult, operator!, result, input, cfg) + + ## update matrix + for j in 1:ndofs4cell_u + ## multiply ansatz function with local jacobian + fill!(tempV, 0) + for d in 1:2 + tempV[1] += jac[1, d] * iduvals[d, j, qp] + tempV[2] += jac[2, d] * iduvals[d, j, qp] + end + for d in 1:4 + tempV[1] += jac[1, 2 + d] * ∇uvals[d, j, qp] + tempV[2] += jac[2, 2 + d] * ∇uvals[d, j, qp] + end + + ## multiply test function operator evaluation + for k in 1:ndofs4cell_u + Aloc[k, j] += dot(tempV, view(iduvals, :, k, qp)) * weights[qp] + end + end + + ## update rhs + mul!(tempV, jac, input) + tempV .-= value + for j in 1:ndofs4cell_u + dof_j = CellDofs_u[j, cell] + b[dof_j] += dot(tempV, view(iduvals, :, j, qp)) * weights[qp] * cellvolumes[cell] + end + end + end + + ## add local matrices to global matrix + Aloc .*= cellvolumes[cell] + for j in 1:ndofs4cell_u + dof_j = CellDofs_u[j, cell] + for k in 1:ndofs4cell_u + dof_k = CellDofs_u[k, cell] + rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k) + end + if (linear) + for k in 1:ndofs4cell_p + dof_k = CellDofs_p[k, cell] + offset_p + rawupdateindex!(A, +, Bloc[j, k], dof_j, dof_k) + rawupdateindex!(A, +, Bloc[j, k], dof_k, dof_j) + end + end + end + fill!(Aloc, 0) + fill!(Bloc, 0) + end + return + end + + function update_system!(linear::Bool, nonlinear::Bool) + barrier(EG, L2G, linear, nonlinear) + return flush!(A) + end + return update_system! end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) - ~, plt = main(; Plotter = Plotter, kwargs...) - scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example210.png"), scene; Plotter = Plotter) + ~, plt = main(; Plotter = Plotter, kwargs...) + scene = GridVisualize.reveal(plt) + return GridVisualize.save(joinpath(dir, "example210.png"), scene; Plotter = Plotter) end end #module diff --git a/examples/Example280_BasisPlotter.jl b/examples/Example280_BasisPlotter.jl index 05132ac1..73d0a9c7 100644 --- a/examples/Example280_BasisPlotter.jl +++ b/examples/Example280_BasisPlotter.jl @@ -18,25 +18,25 @@ using ExtendableGrids ## everything is wrapped in a main function function main(; dim = 1, order = 3) - ## generate two grids - @assert dim in [1, 2] "dim must be 1 or 2" - refgeom = dim == 1 ? Edge1D : Triangle2D - xgrid = reference_domain(refgeom) - - ## set finite element type and get some information - FEType = H1Pk{1, dim, order} - ndofs = get_ndofs(ON_CELLS, FEType, refgeom) - FEType = H1Pk{ndofs, dim, order} - - ## generate FEVector with ncomponents = ndofs - ## that will carry one basis function in each component - FEFunc = FEVector(FESpace{FEType}(xgrid)) - coffsets = ExtendableFEMBase.get_local_coffsets(FEType, ON_CELLS, refgeom) - for j ∈ 1:ndofs - FEFunc[1][j+coffsets[j]] = 1 - end + ## generate two grids + @assert dim in [1, 2] "dim must be 1 or 2" + refgeom = dim == 1 ? Edge1D : Triangle2D + xgrid = reference_domain(refgeom) + + ## set finite element type and get some information + FEType = H1Pk{1, dim, order} + ndofs = get_ndofs(ON_CELLS, FEType, refgeom) + FEType = H1Pk{ndofs, dim, order} + + ## generate FEVector with ncomponents = ndofs + ## that will carry one basis function in each component + FEFunc = FEVector(FESpace{FEType}(xgrid)) + coffsets = ExtendableFEMBase.get_local_coffsets(FEType, ON_CELLS, refgeom) + for j in 1:ndofs + FEFunc[1][j + coffsets[j]] = 1 + end ## plot - println(stdout, unicode_scalarplot(FEFunc[1]; title = "φ", ylim = (-0.5, 1), resolution = dim == 1 ? (40, 10) : (20, 15), nrows = order)) + return println(stdout, unicode_scalarplot(FEFunc[1]; title = "φ", ylim = (-0.5, 1), resolution = dim == 1 ? (40, 10) : (20, 15), nrows = order)) end end diff --git a/examples/Example281_DiscontinuousPlot.jl b/examples/Example281_DiscontinuousPlot.jl index 428445a6..83bd0354 100644 --- a/examples/Example281_DiscontinuousPlot.jl +++ b/examples/Example281_DiscontinuousPlot.jl @@ -19,11 +19,11 @@ using GridVisualize ## function to interpolate function u!(result, qpinfo) - x = qpinfo.x - if qpinfo.region == 1 - result[1] = 2*x[1]*x[2] + x = qpinfo.x + return if qpinfo.region == 1 + result[1] = 2 * x[1] * x[2] elseif qpinfo.region == 2 - result[1] = -1*x[2]*x[1] + 0.5 + result[1] = -1 * x[2] * x[1] + 0.5 else @error "function was evaluated without region information" end @@ -32,21 +32,21 @@ end ## everything is wrapped in a main function function main(; broken = false, nrefs = 3, abs = false, Plotter = nothing) - ## generate two grids - xgrid = grid_unitsquare(Triangle2D) + ## generate two grids + xgrid = grid_unitsquare(Triangle2D) ## mark first two triangles to be in second region xgrid[CellRegions][1:2] .= 2 ## refine - xgrid = uniform_refine(xgrid, nrefs) + xgrid = uniform_refine(xgrid, nrefs) - ## generate coressponding finite element spaces and FEVectors - FES = FESpace{L2P1{1}}(xgrid; broken = broken) - FEFunction = FEVector(FES) + ## generate coressponding finite element spaces and FEVectors + FES = FESpace{L2P1{1}}(xgrid; broken = broken) + FEFunction = FEVector(FES) - ## interpolate function onto first grid - interpolate!(FEFunction[1], u!; bonus_quadorder = 2) + ## interpolate function onto first grid + interpolate!(FEFunction[1], u!; bonus_quadorder = 2) ## get subgrid for each region subgrid1 = subgrid(xgrid, [1]) @@ -60,19 +60,19 @@ function main(; broken = false, nrefs = 3, abs = false, Plotter = nothing) nodevals4nodes1 = nodevalues(FEFunction[1], Identity; abs = abs, regions = [1], nodes = subnodes1) nodevals4nodes2 = nodevalues(FEFunction[1], Identity; abs = abs, regions = [2], nodes = subnodes2) - ## plot + ## plot if Plotter !== nothing - p = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, resolution = (1000, 500)) - gridplot!(p[1,1], xgrid) - scalarplot!(p[1, 2], [subgrid1, subgrid2], xgrid, [view(nodevals4nodes1,:), view(nodevals4nodes2,:)], cellwise = false, levels = 11, title = "u") + p = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, resolution = (1000, 500)) + gridplot!(p[1, 1], xgrid) + scalarplot!(p[1, 2], [subgrid1, subgrid2], xgrid, [view(nodevals4nodes1, :), view(nodevals4nodes2, :)], cellwise = false, levels = 11, title = "u") end return p end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) - plt = main(; Plotter = Plotter, kwargs...) - scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example281.png"), scene; Plotter = Plotter) + plt = main(; Plotter = Plotter, kwargs...) + scene = GridVisualize.reveal(plt) + return GridVisualize.save(joinpath(dir, "example281.png"), scene; Plotter = Plotter) end end diff --git a/examples/Example290_InterpolationBetweenMeshes.jl b/examples/Example290_InterpolationBetweenMeshes.jl index 6946175c..bb3a3e48 100644 --- a/examples/Example290_InterpolationBetweenMeshes.jl +++ b/examples/Example290_InterpolationBetweenMeshes.jl @@ -19,50 +19,50 @@ using GridVisualize ## function to interpolate function u!(result, qpinfo) - x = qpinfo.x - result[1] = sin(4 * pi * x[1]) * sin(4 * pi * x[2]) - result[2] = cos(4 * pi * x[1]) * cos(4 * pi * x[2]) + x = qpinfo.x + result[1] = sin(4 * pi * x[1]) * sin(4 * pi * x[2]) + return result[2] = cos(4 * pi * x[1]) * cos(4 * pi * x[2]) end ## everything is wrapped in a main function -function main(; ν = 1e-3, nrefs = 4, Plotter = nothing) +function main(; ν = 1.0e-3, nrefs = 4, Plotter = nothing) - ## generate two grids - xgrid1 = uniform_refine(grid_unitsquare(Triangle2D), nrefs) - xgrid2 = uniform_refine(xgrid1, 3; store_parents = true) + ## generate two grids + xgrid1 = uniform_refine(grid_unitsquare(Triangle2D), nrefs) + xgrid2 = uniform_refine(xgrid1, 3; store_parents = true) - @show xgrid1 xgrid2 + @show xgrid1 xgrid2 - ## set finite element types for the two grids - FEType1 = H1Pk{2, 2, 2} - FEType2 = H1Pk{2, 2, 1} + ## set finite element types for the two grids + FEType1 = H1Pk{2, 2, 2} + FEType2 = H1Pk{2, 2, 1} - ## generate coressponding finite element spaces and FEVectors - FES1 = FESpace{FEType1}(xgrid1) - FES2 = FESpace{FEType2}(xgrid2) - FEFunction1 = FEVector(FES1) - FEFunction2 = FEVector(FES2) + ## generate coressponding finite element spaces and FEVectors + FES1 = FESpace{FEType1}(xgrid1) + FES2 = FESpace{FEType2}(xgrid2) + FEFunction1 = FEVector(FES1) + FEFunction2 = FEVector(FES2) - ## interpolate function onto first grid - @time interpolate!(FEFunction1[1], u!) - @time interpolate!(FEFunction2[1], u!) + ## interpolate function onto first grid + @time interpolate!(FEFunction1[1], u!) + @time interpolate!(FEFunction2[1], u!) - ## interpolate onto other grid - @time lazy_interpolate!(FEFunction2[1], FEFunction1) - @time lazy_interpolate!(FEFunction2[1], FEFunction1; use_cellparents = true) + ## interpolate onto other grid + @time lazy_interpolate!(FEFunction2[1], FEFunction1) + @time lazy_interpolate!(FEFunction2[1], FEFunction1; use_cellparents = true) - ## plot - p = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (800, 400)) - scalarplot!(p[1, 1], xgrid1, view(nodevalues(FEFunction1[1]), 1, :), levels = 11, title = "u_h ($FEType1, coarse grid)") - scalarplot!(p[1, 2], xgrid2, view(nodevalues(FEFunction2[1]), 1, :), levels = 11, title = "u_h ($FEType2, fine grid)") + ## plot + p = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (800, 400)) + scalarplot!(p[1, 1], xgrid1, view(nodevalues(FEFunction1[1]), 1, :), levels = 11, title = "u_h ($FEType1, coarse grid)") + scalarplot!(p[1, 2], xgrid2, view(nodevalues(FEFunction2[1]), 1, :), levels = 11, title = "u_h ($FEType2, fine grid)") - return p + return p end function generateplots(dir = pwd(); Plotter = nothing, kwargs...) - plt = main(; Plotter = Plotter, kwargs...) - scene = GridVisualize.reveal(plt) - GridVisualize.save(joinpath(dir, "example290.png"), scene; Plotter = Plotter) + plt = main(; Plotter = Plotter, kwargs...) + scene = GridVisualize.reveal(plt) + return GridVisualize.save(joinpath(dir, "example290.png"), scene; Plotter = Plotter) end end diff --git a/pluto-examples/LowLevelNavierStokes.jl b/pluto-examples/LowLevelNavierStokes.jl index 8887f145..0d889420 100644 --- a/pluto-examples/LowLevelNavierStokes.jl +++ b/pluto-examples/LowLevelNavierStokes.jl @@ -6,14 +6,14 @@ using InteractiveUtils # ╔═╡ 014c3380-b361-11ed-273e-37541f1ed74f begin - using ExtendableFEMBase - using ForwardDiff - using DiffResults - using LinearAlgebra - using ExtendableGrids - using ExtendableSparse - using GridVisualize - using PlutoVista + using ExtendableFEMBase + using ForwardDiff + using DiffResults + using LinearAlgebra + using ExtendableGrids + using ExtendableSparse + using GridVisualize + using PlutoVista end # ╔═╡ 3e8fd4f5-b6be-4019-9c76-a80cc985b70e @@ -47,66 +47,66 @@ This is a plot of the computed velocity and pressure: # ╔═╡ 6160364e-6ab2-475d-9345-3436e6f2b3e3 begin - ## PDE data - const μ = 1e-2 - function f!(fval, x, t) # right-hand side - fval[1] = 8.0*π*π*μ * exp(-8.0*π*π*μ*t) * sin(2.0*π*x[1])*sin(2.0*π*x[2]) - fval[2] = 8.0*π*π*μ * exp(-8.0*π*π*μ*t) * cos(2.0*π*x[1])*cos(2.0*π*x[2]) - return nothing - end - - # exact velocity (for boundary data and error calculation) - function u!(uval, qpinfo) - x = qpinfo.x - t = qpinfo.time - uval[1] = exp(-8.0*π*π*μ*t) * sin(2.0*π*x[1])*sin(2.0*π*x[2]) - uval[2] = exp(-8.0*π*π*μ*t) * cos(2.0*π*x[1])*cos(2.0*π*x[2]) - return nothing - end - - ## discretization parameters - const nref = 5 - const teval = 0 - const order = 2 - - ## prepare error calculation - function p!(pval,x,t) # exact pressure (for error calculation) - pval[1] = exp(-16*pi*pi*μ*t)*(cos(4*pi*x[1])-cos(4*pi*x[2]))/4 - return nothing - end + ## PDE data + const μ = 1.0e-2 + function f!(fval, x, t) # right-hand side + fval[1] = 8.0 * π * π * μ * exp(-8.0 * π * π * μ * t) * sin(2.0 * π * x[1]) * sin(2.0 * π * x[2]) + fval[2] = 8.0 * π * π * μ * exp(-8.0 * π * π * μ * t) * cos(2.0 * π * x[1]) * cos(2.0 * π * x[2]) + return nothing + end + + # exact velocity (for boundary data and error calculation) + function u!(uval, qpinfo) + x = qpinfo.x + t = qpinfo.time + uval[1] = exp(-8.0 * π * π * μ * t) * sin(2.0 * π * x[1]) * sin(2.0 * π * x[2]) + uval[2] = exp(-8.0 * π * π * μ * t) * cos(2.0 * π * x[1]) * cos(2.0 * π * x[2]) + return nothing + end + + ## discretization parameters + const nref = 5 + const teval = 0 + const order = 2 + + ## prepare error calculation + function p!(pval, x, t) # exact pressure (for error calculation) + pval[1] = exp(-16 * pi * pi * μ * t) * (cos(4 * pi * x[1]) - cos(4 * pi * x[2])) / 4 + return nothing + end end # ╔═╡ 8e0a11e7-7998-403a-8484-7d0180ea40b2 begin - ## create grid - X = LinRange(0,1,2^nref+1) - Y = LinRange(0,1,2^nref+1) - println("Creating grid...") - @time xgrid = simplexgrid(X,Y) - println("Preparing FaceNodes...") - @time xgrid[FaceNodes] - println("Preparing CellVolumes...") + ## create grid + X = LinRange(0, 1, 2^nref + 1) + Y = LinRange(0, 1, 2^nref + 1) + println("Creating grid...") + @time xgrid = simplexgrid(X, Y) + println("Preparing FaceNodes...") + @time xgrid[FaceNodes] + println("Preparing CellVolumes...") @time xgrid[CellVolumes] - xgrid + xgrid end # ╔═╡ 7e2a082f-5457-47ac-96d5-b688a70a5506 begin - ## create finite element space (Taylor--Hood) - FETypes = [H1Pk{2,2,order}, H1Pk{1,2,order-1}] + ## create finite element space (Taylor--Hood) + FETypes = [H1Pk{2, 2, order}, H1Pk{1, 2, order - 1}] - ## prepare finite element space and dofmaps - println("Creating FESpace...") - @time FES = [FESpace{FETypes[1]}(xgrid; name = "velocity space"), FESpace{FETypes[2]}(xgrid; name = "pressure space")] - FES + ## prepare finite element space and dofmaps + println("Creating FESpace...") + @time FES = [FESpace{FETypes[1]}(xgrid; name = "velocity space"), FESpace{FETypes[2]}(xgrid; name = "pressure space")] + FES end # ╔═╡ 6dcca265-b8fd-4043-a2a7-4ff9bf579dd5 function prepare_assembly!(A, b, FESu, FESp, Solution, f, μ = 1) - A = A.entries - b = b.entries - Solution = Solution.entries + A = A.entries + b = b.entries + Solution = Solution.entries xgrid = FESu.xgrid EG = xgrid[UniqueCellGeometries][1] FEType_u = eltype(FESu) @@ -118,252 +118,254 @@ function prepare_assembly!(A, b, FESu, FESp, Solution, f, μ = 1) ## dofmap CellDofs_u = FESu[ExtendableFEMBase.CellDofs] CellDofs_p = FESp[ExtendableFEMBase.CellDofs] - offset_p = FESu.ndofs - + offset_p = FESu.ndofs + ## quadrature formula - qf = QuadratureRule{Float64, EG}(3*get_polynomialorder(FEType_u, EG)-1) - weights::Vector{Float64} = qf.w - xref::Vector{Vector{Float64}} = qf.xref - nweights::Int = length(weights) - - ## FE basis evaluator - FEBasis_∇u = FEEvaluator(FESu, Gradient, qf) - ∇uvals = FEBasis_∇u.cvals - FEBasis_idu = FEEvaluator(FESu, Identity, qf) - iduvals = FEBasis_idu.cvals - FEBasis_idp = FEEvaluator(FESp, Identity, qf) - idpvals = FEBasis_idp.cvals - - ## prepare automatic differentation of convection operator - function operator!(result, input) - # result = (u ⋅ ∇)u - result[1] = input[1]*input[3]+input[2]*input[4] - result[2] = input[1]*input[5]+input[2]*input[6] - end - result = Vector{Float64}(undef,2) - input = Vector{Float64}(undef,6) + qf = QuadratureRule{Float64, EG}(3 * get_polynomialorder(FEType_u, EG) - 1) + weights::Vector{Float64} = qf.w + xref::Vector{Vector{Float64}} = qf.xref + nweights::Int = length(weights) + + ## FE basis evaluator + FEBasis_∇u = FEEvaluator(FESu, Gradient, qf) + ∇uvals = FEBasis_∇u.cvals + FEBasis_idu = FEEvaluator(FESu, Identity, qf) + iduvals = FEBasis_idu.cvals + FEBasis_idp = FEEvaluator(FESp, Identity, qf) + idpvals = FEBasis_idp.cvals + + ## prepare automatic differentation of convection operator + function operator!(result, input) + # result = (u ⋅ ∇)u + result[1] = input[1] * input[3] + input[2] * input[4] + return result[2] = input[1] * input[5] + input[2] * input[6] + end + result = Vector{Float64}(undef, 2) + input = Vector{Float64}(undef, 6) tempV = zeros(Float64, 2) Dresult = DiffResults.JacobianResult(result, input) cfg = ForwardDiff.JacobianConfig(operator!, result, input, ForwardDiff.Chunk{6}()) jac = DiffResults.jacobian(Dresult) value = DiffResults.value(Dresult) - - + + ## ASSEMBLY LOOP function barrier(EG, L2G::L2GTransformer, linear::Bool, nonlinear::Bool) - ## barrier function to avoid allocations caused by L2G - - ndofs4cell_u::Int = get_ndofs(ON_CELLS, FEType_u, EG) - ndofs4cell_p::Int = get_ndofs(ON_CELLS, FEType_p, EG) - Aloc = zeros(Float64, ndofs4cell_u, ndofs4cell_u) - Bloc = zeros(Float64, ndofs4cell_u, ndofs4cell_p) - dof_j::Int, dof_k::Int = 0, 0 - fval::Vector{Float64} = zeros(Float64,2) + ## barrier function to avoid allocations caused by L2G + + ndofs4cell_u::Int = get_ndofs(ON_CELLS, FEType_u, EG) + ndofs4cell_p::Int = get_ndofs(ON_CELLS, FEType_p, EG) + Aloc = zeros(Float64, ndofs4cell_u, ndofs4cell_u) + Bloc = zeros(Float64, ndofs4cell_u, ndofs4cell_p) + dof_j::Int, dof_k::Int = 0, 0 + fval::Vector{Float64} = zeros(Float64, 2) x::Vector{Float64} = zeros(Float64, 2) - - for cell = 1 : ncells - ## update FE basis evaluators - update_basis!(FEBasis_∇u, cell) - update_basis!(FEBasis_idu, cell) - update_basis!(FEBasis_idp, cell) - - ## assemble local stiffness matrix (symmetric) - if (linear) - for j = 1 : ndofs4cell_u, k = 1 : ndofs4cell_u - temp = 0 - for qp = 1 : nweights - temp += weights[qp] * dot(view(∇uvals,:,j,qp), view(∇uvals,:,k,qp)) - end - Aloc[k,j] = μ * temp - end - - ## assemble div-pressure coupling - for j = 1 : ndofs4cell_u, k = 1 : ndofs4cell_p - temp = 0 - for qp = 1 : nweights - temp -= weights[qp] * (∇uvals[1,j,qp] + ∇uvals[4,j,qp]) * - idpvals[1,k,qp] - end - Bloc[j,k] = temp - end - Bloc .*= cellvolumes[cell] - - ## assemble right-hand side - update_trafo!(L2G, cell) - for j = 1 : ndofs4cell_u - ## right-hand side - temp = 0 - for qp = 1 : nweights - ## get global x for quadrature point - eval_trafo!(x, L2G, xref[qp]) - ## evaluate (f(x), v_j(x)) - f!(fval, x, teval) - temp += weights[qp] * dot(view(iduvals,: , j, qp), fval) - end - ## write into global vector - dof_j = CellDofs_u[j, cell] - b[dof_j] += temp * cellvolumes[cell] - end - end - - ## assemble nonlinear term - if (nonlinear) - for qp = 1 : nweights - fill!(input,0) - for j = 1 : ndofs4cell_u - dof_j = CellDofs_u[j, cell] - for d = 1 : 2 - input[d] += Solution[dof_j] * iduvals[d,j,qp] - end - for d = 1 : 4 - input[2+d] += Solution[dof_j] * ∇uvals[d,j,qp] - end - end - - ## evaluate jacobian - ForwardDiff.chunk_mode_jacobian!(Dresult, operator!, result, input, cfg) - - # update matrix - for j = 1 : ndofs4cell_u - # multiply ansatz function with local jacobian - fill!(tempV,0) - for d = 1 : 2 - tempV[1] += jac[1,d] * iduvals[d,j,qp] - tempV[2] += jac[2,d] * iduvals[d,j,qp] - end - for d = 1 : 4 - tempV[1] += jac[1,2+d] * ∇uvals[d,j,qp] - tempV[2] += jac[2,2+d] * ∇uvals[d,j,qp] - end - - # multiply test function operator evaluation - for k = 1 : ndofs4cell_u - Aloc[k,j] += dot(tempV,view(iduvals,:,k,qp)) * weights[qp] - end - end - - # update rhs - mul!(tempV, jac, input) - tempV .-= value - for j = 1 : ndofs4cell_u - dof_j = CellDofs_u[j, cell] - b[dof_j] += dot(tempV, view(iduvals,:,j,qp)) * weights[qp] * cellvolumes[cell] - end - end - end - - ## add local matrices to global matrix - Aloc .*= cellvolumes[cell] - for j = 1 : ndofs4cell_u - dof_j = CellDofs_u[j, cell] - for k = 1 : ndofs4cell_u - dof_k = CellDofs_u[k, cell] - rawupdateindex!(A, +, Aloc[j,k], dof_j, dof_k) - end - if (linear) - for k = 1 : ndofs4cell_p - dof_k = CellDofs_p[k, cell] + offset_p - rawupdateindex!(A, +, Bloc[j,k], dof_j, dof_k) - rawupdateindex!(A, +, Bloc[j,k], dof_k, dof_j) - end - end - end - fill!(Aloc, 0) - fill!(Bloc, 0) + + for cell in 1:ncells + ## update FE basis evaluators + update_basis!(FEBasis_∇u, cell) + update_basis!(FEBasis_idu, cell) + update_basis!(FEBasis_idp, cell) + + ## assemble local stiffness matrix (symmetric) + if (linear) + for j in 1:ndofs4cell_u, k in 1:ndofs4cell_u + temp = 0 + for qp in 1:nweights + temp += weights[qp] * dot(view(∇uvals, :, j, qp), view(∇uvals, :, k, qp)) + end + Aloc[k, j] = μ * temp + end + + ## assemble div-pressure coupling + for j in 1:ndofs4cell_u, k in 1:ndofs4cell_p + temp = 0 + for qp in 1:nweights + temp -= weights[qp] * (∇uvals[1, j, qp] + ∇uvals[4, j, qp]) * + idpvals[1, k, qp] + end + Bloc[j, k] = temp + end + Bloc .*= cellvolumes[cell] + + ## assemble right-hand side + update_trafo!(L2G, cell) + for j in 1:ndofs4cell_u + ## right-hand side + temp = 0 + for qp in 1:nweights + ## get global x for quadrature point + eval_trafo!(x, L2G, xref[qp]) + ## evaluate (f(x), v_j(x)) + f!(fval, x, teval) + temp += weights[qp] * dot(view(iduvals, :, j, qp), fval) + end + ## write into global vector + dof_j = CellDofs_u[j, cell] + b[dof_j] += temp * cellvolumes[cell] + end + end + + ## assemble nonlinear term + if (nonlinear) + for qp in 1:nweights + fill!(input, 0) + for j in 1:ndofs4cell_u + dof_j = CellDofs_u[j, cell] + for d in 1:2 + input[d] += Solution[dof_j] * iduvals[d, j, qp] + end + for d in 1:4 + input[2 + d] += Solution[dof_j] * ∇uvals[d, j, qp] + end + end + + ## evaluate jacobian + ForwardDiff.chunk_mode_jacobian!(Dresult, operator!, result, input, cfg) + + # update matrix + for j in 1:ndofs4cell_u + # multiply ansatz function with local jacobian + fill!(tempV, 0) + for d in 1:2 + tempV[1] += jac[1, d] * iduvals[d, j, qp] + tempV[2] += jac[2, d] * iduvals[d, j, qp] + end + for d in 1:4 + tempV[1] += jac[1, 2 + d] * ∇uvals[d, j, qp] + tempV[2] += jac[2, 2 + d] * ∇uvals[d, j, qp] + end + + # multiply test function operator evaluation + for k in 1:ndofs4cell_u + Aloc[k, j] += dot(tempV, view(iduvals, :, k, qp)) * weights[qp] + end + end + + # update rhs + mul!(tempV, jac, input) + tempV .-= value + for j in 1:ndofs4cell_u + dof_j = CellDofs_u[j, cell] + b[dof_j] += dot(tempV, view(iduvals, :, j, qp)) * weights[qp] * cellvolumes[cell] + end + end + end + + ## add local matrices to global matrix + Aloc .*= cellvolumes[cell] + for j in 1:ndofs4cell_u + dof_j = CellDofs_u[j, cell] + for k in 1:ndofs4cell_u + dof_k = CellDofs_u[k, cell] + rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k) + end + if (linear) + for k in 1:ndofs4cell_p + dof_k = CellDofs_p[k, cell] + offset_p + rawupdateindex!(A, +, Bloc[j, k], dof_j, dof_k) + rawupdateindex!(A, +, Bloc[j, k], dof_k, dof_j) + end + end + end + fill!(Aloc, 0) + fill!(Bloc, 0) end + return end - function update_system!(linear::Bool, nonlinear::Bool) - barrier(EG, L2G, linear, nonlinear) - flush!(A) - end - update_system! + function update_system!(linear::Bool, nonlinear::Bool) + barrier(EG, L2G, linear, nonlinear) + return flush!(A) + end + return update_system! end # ╔═╡ 7299586b-6859-41af-bcd0-1a0daa454c81 function solve_stokes_lowlevel(FES, μ, f!) - - println("Initializing system...") - Solution = FEVector(FES) - A = FEMatrix(FES) - b = FEVector(FES) - @time update_system! = prepare_assembly!(A, b, FES[1], FES[2], Solution, f!, μ) - @time update_system!(true, false) - Alin = deepcopy(A) # = keep linear part of system matrix - blin = deepcopy(b) # = keep linear part of right-hand side - - println("Prepare boundary conditions...") - @time begin - u_init = FEVector(FES) - interpolate!(u_init[1], u!; time = teval) - - fixed_dofs = [size(A.entries,1)] # fix one pressure dof = last dof - BFaceDofs::Adjacency{Int32} = FES[1][ExtendableFEMBase.BFaceDofs] - nbfaces::Int = num_sources(BFaceDofs) - AM::ExtendableSparseMatrix{Float64,Int64} = A.entries - dof_j::Int = 0 - for bface = 1 : nbfaces - for j = 1 : num_targets(BFaceDofs,1) - dof_j = BFaceDofs[j, bface] - push!(fixed_dofs, dof_j) - end - end - end - - - for it = 1 : 20 - ## solve - println("\nITERATION $it\n=============") - println("Solving linear system...") - @time copyto!(Solution.entries, A.entries \ b.entries) - res = A.entries.cscmatrix * Solution.entries .- b.entries - for dof in fixed_dofs - res[dof] = 0 - end - linres = norm(res) - println("linear residual = $linres") - - fill!(A.entries.cscmatrix.nzval,0) - fill!(b.entries,0) - println("Updating linear system...") - @time begin - update_system!(false,true) - A.entries.cscmatrix += Alin.entries.cscmatrix - b.entries .+= blin.entries - end - - ## fix boundary dofs - for dof in fixed_dofs - AM[dof,dof] = 1e60 - b.entries[dof] = 1e60 * u_init.entries[dof] - end - ExtendableSparse.flush!(A.entries) - - ## calculate nonlinear residual - res = A.entries.cscmatrix * Solution.entries .- b.entries - for dof in fixed_dofs - res[dof] = 0 - end - nlres = norm(res) - println("nonlinear residual = $nlres") - if nlres < max(1e-12, 20*linres) - break - end - end - - return Solution, u_init + + println("Initializing system...") + Solution = FEVector(FES) + A = FEMatrix(FES) + b = FEVector(FES) + @time update_system! = prepare_assembly!(A, b, FES[1], FES[2], Solution, f!, μ) + @time update_system!(true, false) + Alin = deepcopy(A) # = keep linear part of system matrix + blin = deepcopy(b) # = keep linear part of right-hand side + + println("Prepare boundary conditions...") + @time begin + u_init = FEVector(FES) + interpolate!(u_init[1], u!; time = teval) + + fixed_dofs = [size(A.entries, 1)] # fix one pressure dof = last dof + BFaceDofs::Adjacency{Int32} = FES[1][ExtendableFEMBase.BFaceDofs] + nbfaces::Int = num_sources(BFaceDofs) + AM::ExtendableSparseMatrix{Float64, Int64} = A.entries + dof_j::Int = 0 + for bface in 1:nbfaces + for j in 1:num_targets(BFaceDofs, 1) + dof_j = BFaceDofs[j, bface] + push!(fixed_dofs, dof_j) + end + end + end + + + for it in 1:20 + ## solve + println("\nITERATION $it\n=============") + println("Solving linear system...") + @time copyto!(Solution.entries, A.entries \ b.entries) + res = A.entries.cscmatrix * Solution.entries .- b.entries + for dof in fixed_dofs + res[dof] = 0 + end + linres = norm(res) + println("linear residual = $linres") + + fill!(A.entries.cscmatrix.nzval, 0) + fill!(b.entries, 0) + println("Updating linear system...") + @time begin + update_system!(false, true) + A.entries.cscmatrix += Alin.entries.cscmatrix + b.entries .+= blin.entries + end + + ## fix boundary dofs + for dof in fixed_dofs + AM[dof, dof] = 1.0e60 + b.entries[dof] = 1.0e60 * u_init.entries[dof] + end + ExtendableSparse.flush!(A.entries) + + ## calculate nonlinear residual + res = A.entries.cscmatrix * Solution.entries .- b.entries + for dof in fixed_dofs + res[dof] = 0 + end + nlres = norm(res) + println("nonlinear residual = $nlres") + if nlres < max(1.0e-12, 20 * linres) + break + end + end + + return Solution, u_init end # ╔═╡ 0785624c-e95a-4e76-8071-2eea591091e0 begin - ## call low level solver - sol, u_init = solve_stokes_lowlevel(FES, μ, f!) - sol + ## call low level solver + sol, u_init = solve_stokes_lowlevel(FES, μ, f!) + sol end # ╔═╡ 6f7a1407-dfb3-497b-8c57-8efac6592194 -[tricontour(xgrid[Coordinates],xgrid[CellNodes],nodevalues(sol[1]; abs = true)[:]; levels = 5), -tricontour(xgrid[Coordinates],xgrid[CellNodes],nodevalues(sol[2])[:]; levels = 5) +[ + tricontour(xgrid[Coordinates], xgrid[CellNodes], nodevalues(sol[1]; abs = true)[:]; levels = 5), + tricontour(xgrid[Coordinates], xgrid[CellNodes], nodevalues(sol[2])[:]; levels = 5), ] diff --git a/pluto-examples/LowLevelPoisson.jl b/pluto-examples/LowLevelPoisson.jl index 0210da1b..9cd6000c 100644 --- a/pluto-examples/LowLevelPoisson.jl +++ b/pluto-examples/LowLevelPoisson.jl @@ -6,11 +6,11 @@ using InteractiveUtils # ╔═╡ 014c3380-b361-11ed-273e-37541f1ed74f begin - using ExtendableFEMBase - using ExtendableGrids - using ExtendableSparse - using GridVisualize - using PlutoVista + using ExtendableFEMBase + using ExtendableGrids + using ExtendableSparse + using GridVisualize + using PlutoVista GridVisualize.default_plotter!(PlutoVista) end @@ -37,38 +37,38 @@ The weak formulation seeks ``u \in V := H^1_0(\Omega)`` such that # ╔═╡ 6160364e-6ab2-475d-9345-3436e6f2b3e3 begin - ## PDE data - μ = 1.0 - f = x -> x[1] - x[2] + ## PDE data + μ = 1.0 + f = x -> x[1] - x[2] - ## discretization parameters - nref = 9 - order = 1 + ## discretization parameters + nref = 9 + order = 1 end # ╔═╡ 8e0a11e7-7998-403a-8484-7d0180ea40b2 begin - ## create grid - X = LinRange(0,1,2^nref+1) - Y = LinRange(0,1,2^nref+1) - println("Creating grid...") - @time xgrid = simplexgrid(X,Y) - println("Preparing FaceNodes...") - @time xgrid[FaceNodes] - println("Preparing CellVolumes...") + ## create grid + X = LinRange(0, 1, 2^nref + 1) + Y = LinRange(0, 1, 2^nref + 1) + println("Creating grid...") + @time xgrid = simplexgrid(X, Y) + println("Preparing FaceNodes...") + @time xgrid[FaceNodes] + println("Preparing CellVolumes...") @time xgrid[CellVolumes] - xgrid + xgrid end # ╔═╡ 7e2a082f-5457-47ac-96d5-b688a70a5506 begin - ## create finite element space - FEType = H1Pk{1,2,order} + ## create finite element space + FEType = H1Pk{1, 2, order} - ## prepare finite element space and dofmaps - println("Creating FESpace...") + ## prepare finite element space and dofmaps + println("Creating FESpace...") @time FES = FESpace{FEType}(xgrid) - FES + FES end # ╔═╡ 6dcca265-b8fd-4043-a2a7-4ff9bf579dd5 @@ -81,126 +81,127 @@ function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1) ## dofmap CellDofs = FES[ExtendableFEMBase.CellDofs] - + ## quadrature formula - qf = QuadratureRule{Float64, EG}(2*(get_polynomialorder(FEType, EG)-1)) - weights::Vector{Float64} = qf.w - xref::Vector{Vector{Float64}} = qf.xref - nweights::Int = length(weights) - - ## FE basis evaluator - FEBasis_∇ = FEEvaluator(FES, Gradient, qf) - ∇vals = FEBasis_∇.cvals - FEBasis_id = FEEvaluator(FES, Identity, qf) - idvals = FEBasis_id.cvals - - + qf = QuadratureRule{Float64, EG}(2 * (get_polynomialorder(FEType, EG) - 1)) + weights::Vector{Float64} = qf.w + xref::Vector{Vector{Float64}} = qf.xref + nweights::Int = length(weights) + + ## FE basis evaluator + FEBasis_∇ = FEEvaluator(FES, Gradient, qf) + ∇vals = FEBasis_∇.cvals + FEBasis_id = FEEvaluator(FES, Identity, qf) + idvals = FEBasis_id.cvals + + cellvolumes = xgrid[CellVolumes] - + ## ASSEMBLY LOOP function barrier(EG, L2G::L2GTransformer) - ## barrier function to avoid allocations by EG dispatch - - ndofs4cell::Int = get_ndofs(ON_CELLS, FEType, EG) - Aloc = zeros(Float64, ndofs4cell, ndofs4cell) + ## barrier function to avoid allocations by EG dispatch + + ndofs4cell::Int = get_ndofs(ON_CELLS, FEType, EG) + Aloc = zeros(Float64, ndofs4cell, ndofs4cell) ncells::Int = num_cells(xgrid) - dof_j::Int, dof_k::Int = 0, 0 + dof_j::Int, dof_k::Int = 0, 0 x::Vector{Float64} = zeros(Float64, 2) - - for cell = 1 : ncells - ## update FE basis evaluators - FEBasis_∇.citem[] = cell - update_basis!(FEBasis_∇) - - ## assemble local stiffness matrix - for j = 1 : ndofs4cell, k = j : ndofs4cell - temp = 0 - for qp = 1 : nweights - temp += weights[qp] * dot(view(∇vals,:,j,qp), view(∇vals,:,k,qp)) - end - Aloc[j,k] = temp - end - Aloc .*= μ * cellvolumes[cell] - - ## add local matrix to global matrix - for j = 1 : ndofs4cell - dof_j = CellDofs[j, cell] - for k = j : ndofs4cell - dof_k = CellDofs[k, cell] - if abs(Aloc[j,k]) > 1e-15 - # write into sparse matrix, only lines with allocations - rawupdateindex!(A, +, Aloc[j,k], dof_j, dof_k) - if k > j - rawupdateindex!(A, +, Aloc[j,k], dof_k, dof_j) - end - end - end - end - fill!(Aloc, 0) - - ## assemble right-hand side + + for cell in 1:ncells + ## update FE basis evaluators + FEBasis_∇.citem[] = cell + update_basis!(FEBasis_∇) + + ## assemble local stiffness matrix + for j in 1:ndofs4cell, k in j:ndofs4cell + temp = 0 + for qp in 1:nweights + temp += weights[qp] * dot(view(∇vals, :, j, qp), view(∇vals, :, k, qp)) + end + Aloc[j, k] = temp + end + Aloc .*= μ * cellvolumes[cell] + + ## add local matrix to global matrix + for j in 1:ndofs4cell + dof_j = CellDofs[j, cell] + for k in j:ndofs4cell + dof_k = CellDofs[k, cell] + if abs(Aloc[j, k]) > 1.0e-15 + # write into sparse matrix, only lines with allocations + rawupdateindex!(A, +, Aloc[j, k], dof_j, dof_k) + if k > j + rawupdateindex!(A, +, Aloc[j, k], dof_k, dof_j) + end + end + end + end + fill!(Aloc, 0) + + ## assemble right-hand side update_trafo!(L2G, cell) - for j = 1 : ndofs4cell + for j in 1:ndofs4cell ## right-hand side temp = 0 - for qp = 1 : nweights + for qp in 1:nweights ## get global x for quadrature point eval_trafo!(x, L2G, xref[qp]) ## evaluate (f(x), v_j(x)) temp += weights[qp] * idvals[1, j, qp] * f(x) end - ## write into global vector + ## write into global vector dof_j = CellDofs[j, cell] b[dof_j] += temp * cellvolumes[cell] end end + return end barrier(EG, L2G) - flush!(A) + return flush!(A) end # ╔═╡ 7299586b-6859-41af-bcd0-1a0daa454c81 function solve_poisson_lowlevel(FES, μ, f) - - Solution = FEVector(FES) - FES = Solution[1].FES - A = FEMatrix(FES, FES) - b = FEVector(FES) - println("Assembling operators...") - @time assemble!(A.entries, b.entries, FES, f, μ) + + Solution = FEVector(FES) + FES = Solution[1].FES + A = FEMatrix(FES, FES) + b = FEVector(FES) + println("Assembling operators...") + @time assemble!(A.entries, b.entries, FES, f, μ) ## fix boundary dofs - println("Assembling boundary data...") - @time begin - BFaceDofs::Adjacency{Int32} = FES[ExtendableFEMBase.BFaceDofs] - nbfaces::Int = num_sources(BFaceDofs) - AM::ExtendableSparseMatrix{Float64,Int64} = A.entries - dof_j::Int = 0 - for bface = 1 : nbfaces - for j = 1 : num_targets(BFaceDofs,1) - dof_j = BFaceDofs[j, bface] - AM[dof_j,dof_j] = 1e60 - b.entries[dof_j] = 0 - end - end - end - ExtendableSparse.flush!(A.entries) + println("Assembling boundary data...") + @time begin + BFaceDofs::Adjacency{Int32} = FES[ExtendableFEMBase.BFaceDofs] + nbfaces::Int = num_sources(BFaceDofs) + AM::ExtendableSparseMatrix{Float64, Int64} = A.entries + dof_j::Int = 0 + for bface in 1:nbfaces + for j in 1:num_targets(BFaceDofs, 1) + dof_j = BFaceDofs[j, bface] + AM[dof_j, dof_j] = 1.0e60 + b.entries[dof_j] = 0 + end + end + end + ExtendableSparse.flush!(A.entries) ## solve - println("Solving linear system...") + println("Solving linear system...") @time copyto!(Solution.entries, A.entries \ b.entries) - return Solution + return Solution end # ╔═╡ 0785624c-e95a-4e76-8071-2eea591091e0 begin - ## call low level solver - sol = solve_poisson_lowlevel(FES, μ, f) + ## call low level solver + sol = solve_poisson_lowlevel(FES, μ, f) end # ╔═╡ 6f7a1407-dfb3-497b-8c57-8efac6592194 -tricontour(xgrid[Coordinates],xgrid[CellNodes],sol.entries[1:num_nodes(xgrid)]; levels = 5, resolution = (500,500)) +tricontour(xgrid[Coordinates], xgrid[CellNodes], sol.entries[1:num_nodes(xgrid)]; levels = 5, resolution = (500, 500)) # ╔═╡ 00000000-0000-0000-0000-000000000001 PLUTO_PROJECT_TOML_CONTENTS = """ diff --git a/src/fedefs/h1_bubble.jl b/src/fedefs/h1_bubble.jl index af17b47c..af7e159d 100644 --- a/src/fedefs/h1_bubble.jl +++ b/src/fedefs/h1_bubble.jl @@ -38,7 +38,7 @@ interior_dofs_offset(::Type{ON_CELLS}, ::Type{H1BUBBLE{ncomponents}}, EG::Type{< init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: H1BUBBLE, APT} = MomentInterpolator(FES, ON_CELLS) function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: H1BUBBLE, APT} - get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items, kwargs...) + return get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items, kwargs...) end function get_basis(::Type{<:AssemblyType}, ::Type{H1BUBBLE{ncomponents}}, ::Type{<:AbstractElementGeometry1D}) where {ncomponents} diff --git a/src/fedefs/h1_p2b.jl b/src/fedefs/h1_p2b.jl index 02f461a0..17ce2461 100644 --- a/src/fedefs/h1_p2b.jl +++ b/src/fedefs/h1_p2b.jl @@ -95,7 +95,7 @@ function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, end # fix cell bubble value by preserving integral mean - get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...) + return get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...) end function get_basis(AT::Union{Type{<:ON_FACES}, Type{<:ON_BFACES}}, ::Type{H1P2B{ncomponents, edim}}, EG::Type{<:AbstractElementGeometry}) where {ncomponents, edim} diff --git a/src/fedefs/h1v_p1teb.jl b/src/fedefs/h1v_p1teb.jl index 3b37a201..e96e7ef1 100644 --- a/src/fedefs/h1v_p1teb.jl +++ b/src/fedefs/h1v_p1teb.jl @@ -49,7 +49,7 @@ end function TEB_tangentflux_eval_3d!(grid) edgetangents = grid[EdgeTangents] function closure(result, f, qpinfo) - result[1] = dot(f, view(edgetangents, :, qpinfo.item)) + return result[1] = dot(f, view(edgetangents, :, qpinfo.item)) end return closure end diff --git a/src/fedefs/hcurl_n0.jl b/src/fedefs/hcurl_n0.jl index fb384b4d..12773428 100644 --- a/src/fedefs/hcurl_n0.jl +++ b/src/fedefs/hcurl_n0.jl @@ -47,7 +47,7 @@ init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {T function N0_tangentflux_eval_3d!(grid) edgetangents = grid[EdgeTangents] function closure(result, f, qpinfo) - result[1] = dot(f, view(edgetangents, :, qpinfo.item)) + return result[1] = dot(f, view(edgetangents, :, qpinfo.item)) end return closure end diff --git a/src/fedefs/hdiv_bdm1.jl b/src/fedefs/hdiv_bdm1.jl index 71f8f93d..87ea2dc6 100644 --- a/src/fedefs/hdiv_bdm1.jl +++ b/src/fedefs/hdiv_bdm1.jl @@ -41,10 +41,10 @@ isdefined(FEType::Type{<:HDIVBDM1}, ::Type{<:Quadrilateral2D}) = true isdefined(FEType::Type{<:HDIVBDM1}, ::Type{<:Tetrahedron3D}) = true function BDM1_normalflux_eval!(dim) - function closure(result, f, qpinfo) + return function closure(result, f, qpinfo) result[1] = dot(f, qpinfo.normal) result[2] = result[1] * (qpinfo.xref[1] - 1 // dim) - if dim == 3 + return if dim == 3 result[3] = result[1] * (qpinfo.xref[2] - 1 // dim) end end @@ -53,7 +53,7 @@ init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {T function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM1, APT} - get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) + return get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) end function ExtendableGrids.interpolate!(Target, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, data; items = [], kwargs...) where {Tv, Ti, FEType <: HDIVBDM1, APT} diff --git a/src/fedefs/hdiv_bdm2.jl b/src/fedefs/hdiv_bdm2.jl index 0275ac05..1f2c6647 100644 --- a/src/fedefs/hdiv_bdm2.jl +++ b/src/fedefs/hdiv_bdm2.jl @@ -32,18 +32,18 @@ interior_dofs_offset(::Type{<:ON_CELLS}, ::Type{<:HDIVBDM2{2}}, ::Type{<:Triangl function BDM2_normalflux_eval!(dim) @assert dim == 2 "BDM3 for dim=3 not available yet" - function closure(result, f, qpinfo) + return function closure(result, f, qpinfo) result[1] = dot(f, qpinfo.normal) result[2] = result[1] * (qpinfo.xref[1] - 1 // dim) - result[3] = result[1] * (qpinfo.xref[1]^2 - qpinfo.xref[1] + 1 // 6) + return result[3] = result[1] * (qpinfo.xref[1]^2 - qpinfo.xref[1] + 1 // 6) end end init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HDIVBDM2, APT} = FunctionalInterpolator(BDM2_normalflux_eval!(FEType.parameters[1]), FES, ON_FACES; bonus_quadorder = 2) -init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: HDIVBDM2{2}, APT} = MomentInterpolator(FES, ON_CELLS; FEType_moments = H1MINI{1,2}, moments_dofs = [1,2,4], moments_operator = CurlScalar) +init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: HDIVBDM2{2}, APT} = MomentInterpolator(FES, ON_CELLS; FEType_moments = H1MINI{1, 2}, moments_dofs = [1, 2, 4], moments_operator = CurlScalar) function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM2, APT} - get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) + return get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) end function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVBDM2, APT} @@ -53,7 +53,7 @@ function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{T # set interior dofs such that moments with ∇P1 and curl(b) are preserved # (b is the cell bubble) - get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...) + return get_interpolator(FE, ON_CELLS).evaluate!(Target, exact_function!, items; kwargs...) end ## only normalfluxes on faces diff --git a/src/fedefs/hdiv_rt1.jl b/src/fedefs/hdiv_rt1.jl index 66b73067..156dcea4 100644 --- a/src/fedefs/hdiv_rt1.jl +++ b/src/fedefs/hdiv_rt1.jl @@ -41,10 +41,10 @@ interior_dofs_offset(::Type{<:ON_CELLS}, ::Type{<:HDIVRT1{2}}, ::Type{<:Triangle interior_dofs_offset(::Type{<:ON_CELLS}, ::Type{<:HDIVRT1{3}}, ::Type{<:Tetrahedron3D}) = 12 function RT1_normalflux_eval!(dim) - function closure(result, f, qpinfo) + return function closure(result, f, qpinfo) result[1] = dot(f, qpinfo.normal) result[2] = result[1] * (qpinfo.xref[1] - 1 // dim) - if dim == 3 + return if dim == 3 result[3] = result[1] * (qpinfo.xref[2] - 1 // dim) end end @@ -53,7 +53,7 @@ init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {T init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {Tv, Ti, FEType <: HDIVRT1, APT} = MomentInterpolator(FES, ON_CELLS) function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVRT1, APT} - get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) + return get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) end function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, FEType <: HDIVRT1, APT} diff --git a/src/fedefs/hdiv_rtk.jl b/src/fedefs/hdiv_rtk.jl index 1fb9c0d0..c84a894c 100644 --- a/src/fedefs/hdiv_rtk.jl +++ b/src/fedefs/hdiv_rtk.jl @@ -33,11 +33,12 @@ interior_dofs_offset(::Type{<:ON_CELLS}, ::Type{<:HDIVRTk{2, order}}, ::Type{<:T function RTk_normalflux_eval!(order) moments_weights = basis.(ShiftedLegendre, (0:order)) - function closure(result, f, qpinfo) + return function closure(result, f, qpinfo) result[1] = dot(f, qpinfo.normal) for j in 1:order result[j + 1] = result[1] * moments_weights[j + 1](qpinfo.xref[1]) end + return end end init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_FACES}) where {Tv, Ti, FEType <: HDIVRTk, APT} = FunctionalInterpolator(RTk_normalflux_eval!(FEType.parameters[2]), FES, ON_FACES; bonus_quadorder = FEType.parameters[2]) @@ -45,7 +46,7 @@ init_interpolator!(FES::FESpace{Tv, Ti, FEType, APT}, ::Type{ON_CELLS}) where {T function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, HDIVRTk{edim, order}, APT}, ::Type{ON_FACES}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, edim, order, APT} - get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) + return get_interpolator(FE, ON_FACES).evaluate!(Target, exact_function!, items; kwargs...) end function ExtendableGrids.interpolate!(Target::AbstractArray{T, 1}, FE::FESpace{Tv, Ti, HDIVRTk{edim, order}, APT}, ::Type{ON_CELLS}, exact_function!; items = [], kwargs...) where {T, Tv, Ti, edim, order, APT} diff --git a/src/fematrix.jl b/src/fematrix.jl index 6ce4d467..ac9993f5 100644 --- a/src/fematrix.jl +++ b/src/fematrix.jl @@ -32,7 +32,7 @@ $(TYPEDSIGNATURES) Custom `show` function for `FEMatrixBlock` that prints its coordinates and the name. """ function Base.show(io::IO, FEB::FEMatrixBlock) - @printf(io, "[%d:%d,%d:%d]: %s", FEB.offset+1, FEB.last_index, FEB.offsetY+1, FEB.last_indexY, FEB.name) + @printf(io, "[%d:%d,%d:%d]: %s", FEB.offset + 1, FEB.last_index, FEB.offsetY + 1, FEB.last_indexY, FEB.name) return nothing end @@ -93,11 +93,10 @@ Base.getindex(FEF::FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtotal}, i::Int, Base.getindex(FEB::FEMatrixBlock, i::Int, j::Int) = FEB.entries[FEB.offset + i, FEB.offsetY + j] Base.getindex(FEB::FEMatrixBlock, i::Any, j::Any) = FEB.entries[FEB.offset .+ i, FEB.offsetY .+ j] Base.setindex!(FEB::FEMatrixBlock, v, i::Int, j::Int) = setindex!(FEB.entries, v, FEB.offset + i, FEB.offsetY + j) -Base.first(FEB::FEMatrixBlock) = (FEB.offset+1, FEB.offsetY+1) +Base.first(FEB::FEMatrixBlock) = (FEB.offset + 1, FEB.offsetY + 1) Base.last(FEB::FEMatrixBlock) = (FEB.last_index, FEB.last_indexY) - """ $(TYPEDSIGNATURES) @@ -148,7 +147,7 @@ function Base.show(io::IO, FEM::FEMatrix{TvM, TiM, TvG, TiG, nbrow, nbcol, nbtot n = mod(j - 1, nbrow) + 1 m = Int(ceil(j / nbrow)) @printf(io, " [%2d,%2d] |", m, n) - @printf(io, " (%6d,%6d) |", FEM[j].offset+1, FEM[j].offsetY+1) + @printf(io, " (%6d,%6d) |", FEM[j].offset + 1, FEM[j].offsetY + 1) @printf(io, " (%6d,%6d) |", FEM[j].last_index, FEM[j].last_indexY) @printf(io, " (%6d,%6d) |", FEM[j].FES.ndofs, FEM[j].FESY.ndofs) @printf(io, " %s \n", FEM[j].name) @@ -624,8 +623,8 @@ function submatrix(A::AbstractExtendableSparseMatrixCSC{Tv, Ti}, srows, scols; f cscmat::SparseMatrixCSC{Tv, Ti} = A.cscmatrix rows::Array{Ti, 1} = rowvals(cscmat) valsA = cscmat.nzval - nrowsA = size(A,1) - ncolsA = size(A,2) + nrowsA = size(A, 1) + ncolsA = size(A, 2) S = ExtendableSparseMatrix{Tv, Ti}(length(srows), length(scols)) @assert maximum(srows) <= size(A, 1) "rows exceeds rowcount of A" @assert maximum(scols) <= size(A, 2) "cols exceeds colcount of A" @@ -641,7 +640,7 @@ function submatrix(A::AbstractExtendableSparseMatrixCSC{Tv, Ti}, srows, scols; f for r in nzrange(cscmat, scol) if newrows[rows[r]] > 0 _addnz(S, newrows[rows[r]], newcols[scol], valsA[r], factor) - else + else break end end @@ -656,5 +655,5 @@ $(TYPEDSIGNATURES) Returns the FEMatrixBlock as an ExtendableSparseMatrix """ function submatrix(A::FEMatrixBlock{Tv, Ti}) where {Tv, Ti} - return submatrix(A.entries, A.offset+1:A.last_index, A.offsetY+1:A.last_indexY) + return submatrix(A.entries, (A.offset + 1):A.last_index, (A.offsetY + 1):A.last_indexY) end diff --git a/src/fevector.jl b/src/fevector.jl index cfcd314a..bb51b300 100644 --- a/src/fevector.jl +++ b/src/fevector.jl @@ -35,8 +35,8 @@ $(TYPEDSIGNATURES) Custom `show` function for `FEVectorBlock` that prints some information and the view of that block. """ function Base.show(io::IO, ::MIME"text/plain", FEB::FEVectorBlock) - @printf(io, "block %s [%d:%d] = ", FEB.name, FEB.offset+1, FEB.last_index) - show(io, view(FEB)) + @printf(io, "block %s [%d:%d] = ", FEB.name, FEB.offset + 1, FEB.last_index) + return show(io, view(FEB)) end """ @@ -69,7 +69,7 @@ Base.setindex!(FEB::FEVectorBlock, v, i::AbstractArray) = (FEB.entries[FEB.offse Base.eltype(::FEVector{T}) where {T} = T Base.size(FEF::FEVector) = size(FEF.FEVectorBlocks) Base.size(FEB::FEVectorBlock) = FEB.last_index - FEB.offset -Base.first(FEB::FEVectorBlock) = FEB.offset+1 +Base.first(FEB::FEVectorBlock) = FEB.offset + 1 Base.last(FEB::FEVectorBlock) = FEB.last_index @@ -197,7 +197,7 @@ function Base.show(io::IO, FEF::FEVector) print(io, " block | starts | ends | length | min / max \t| FEType \t\t (name/tag)") for j in 1:length(FEF) @printf(io, "\n [%5d] |", j) - @printf(io, " %6d |", FEF[j].offset+1) + @printf(io, " %6d |", FEF[j].offset + 1) @printf(io, " %6d |", FEF[j].last_index) @printf(io, " %6d |", FEF[j].FES.ndofs) ext = extrema(view(FEF[j])) diff --git a/src/finiteelements.jl b/src/finiteelements.jl index e51469b6..b616ac76 100644 --- a/src/finiteelements.jl +++ b/src/finiteelements.jl @@ -490,7 +490,7 @@ function get_reconstruction_matrix(T::Type{<:Real}, FE::FESpace, FER::FESpace) for j in 2:length(EG) append!(rhandlers, [get_reconstruction_coefficients(ON_CELLS, FE, FER, EG[j])]) append!(chandlers, [get_coefficients(ON_CELLS, FER, EG[j], xgrid)]) - append!(shandlers, [get_basissubset(ON_CELLS, FER, EG[j],)]) + append!(shandlers, [get_basissubset(ON_CELLS, FER, EG[j])]) end ndofs_FE = zeros(Int, length(EG)) diff --git a/src/interpolators.jl b/src/interpolators.jl index e4c7c07e..211d6b9b 100644 --- a/src/interpolators.jl +++ b/src/interpolators.jl @@ -27,11 +27,12 @@ function at the requested nodes (items) of the grid. In broken mode only ON_CELL refer to cells and all nodes of these cells are evaluated. """ function NodalInterpolator( - FES::FESpace{T}, - grid = FES.dofgrid; - broken = FES.broken, - component_offset = FES.coffset, - kwargs...) where {T} + FES::FESpace{T}, + grid = FES.dofgrid; + broken = FES.broken, + component_offset = FES.coffset, + kwargs... + ) where {T} FEType = eltype(FES) ncomponents = get_ncomponents(FEType) @@ -49,7 +50,7 @@ function NodalInterpolator( QP.time = time QP.params = params === nothing ? [] : params if isempty(items) - items = 1 : ncells + items = 1:ncells end for cell in items if cell < 1 @@ -68,6 +69,7 @@ function NodalInterpolator( end end end + return end return NodalInterpolator(evaluate_broken!) else @@ -79,7 +81,7 @@ function NodalInterpolator( QP.time = time QP.params = params === nothing ? [] : params if isempty(items) - items = 1 : nnodes + items = 1:nnodes end for j in items if j < 1 @@ -95,6 +97,7 @@ function NodalInterpolator( target[j + offset4component[k]] = result[k] end end + return end return NodalInterpolator(evaluate!) end @@ -117,19 +120,20 @@ basis functions of the moments of type FEType_ref and with moments_operator. In the mass matrix instead is formed from the scalar product of the interior basis functions. """ function MomentInterpolator( - FE::FESpace{Tv, Ti, FEType, APT}, - AT::Type{<:AssemblyType}, - xgrid = FE.dofgrid; - operator = Identity, - FEType_ref = :auto, - FEType_moments = :auto, - moments_operator = operator, - moments_dofs = Int[], - bestapprox = false, - order = 0, - coffset::Int = -1, - componentwise = true, - kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT} + FE::FESpace{Tv, Ti, FEType, APT}, + AT::Type{<:AssemblyType}, + xgrid = FE.dofgrid; + operator = Identity, + FEType_ref = :auto, + FEType_moments = :auto, + moments_operator = operator, + moments_dofs = Int[], + bestapprox = false, + order = 0, + coffset::Int = -1, + componentwise = true, + kwargs... + ) where {Tv, Ti, FEType <: AbstractFiniteElement, APT} itemvolumes = xgrid[GridComponentVolumes4AssemblyType(AT)] itemnodes = xgrid[GridComponentNodes4AssemblyType(AT)] @@ -189,9 +193,9 @@ function MomentInterpolator( ## check if mass matrix can be computed once or needs to be recomputed on every mesh if FEType_ref <: AbstractH1FiniteElement && - !(FEType_ref <: AbstractH1FiniteElementWithCoefficients) && - FEType_moments <: AbstractH1FiniteElement && - moments_operator == Identity + !(FEType_ref <: AbstractH1FiniteElementWithCoefficients) && + FEType_moments <: AbstractH1FiniteElement && + moments_operator == Identity fixed_mass_matrix = true else fixed_mass_matrix = false @@ -200,15 +204,15 @@ function MomentInterpolator( moments_basis! = get_basis(AT, FEType_moments, EG) nmoments::Int = get_ndofs_all(AT, FEType_moments, EG) if isempty(moments_dofs) - moments_dofs = Array{Int,1}(1:nmoments) + moments_dofs = Array{Int, 1}(1:nmoments) else nmoments = length(moments_dofs) end xgrid_ref = reference_domain(EG) idofs = zeros(Int, 0) - if bestapprox + if bestapprox FEType_moments = FEType_ref - append!(idofs, interior_offset+1:get_ndofs(AT, FEType_ref, EG)) + append!(idofs, (interior_offset + 1):get_ndofs(AT, FEType_ref, EG)) nmoments = length(idofs) else if FEType_ref <: AbstractH1FiniteElement && componentwise @@ -303,7 +307,7 @@ function MomentInterpolator( # prepare integration of moments L2G = L2GTransformer(EG, xgrid, AT) - current_quadorder = 2*order_FE + current_quadorder = 2 * order_FE QF = QuadratureRule{Tv, EG}(current_quadorder) f_moments = zeros(Tv, nmoments) result_f = zeros(Tv, ncomponents) @@ -337,7 +341,7 @@ function MomentInterpolator( if moments_operator !== Identity update_basis!(FEB_moments, item) end - + ## integrate moments of function fill!(f_moments, 0) for qp in 1:nweights @@ -380,7 +384,7 @@ function MomentInterpolator( for dof in 1:interior_offset for i in 1:nweights eval_febe!(basisval, FEB, dof, i) - for m in 1:nmoments, k = 1:ncomponents + for m in 1:nmoments, k in 1:ncomponents f_moments[m] -= basisval[k] * FEB_moments.cvals[k, moments_dofs[m], i] * target[itemdofs[dof, item]] * weights[i] end end @@ -456,17 +460,18 @@ The functionals are corrected by the operator evaluations of the fixed dofs. The mean parameter decides if the functional is divided by the area in the end (if mean = true). """ function FunctionalInterpolator( - functionals!::Function, - FE::FESpace{Tv, Ti, FEType, APT}, - AT::Type{<:AssemblyType} = ON_FACES, - xgrid = FE.dofgrid; - bonus_quadorder = 0, - operator = NormalFlux, - nfluxes = 0, - dofs = [], - mean = false, - kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT} - + functionals!::Function, + FE::FESpace{Tv, Ti, FEType, APT}, + AT::Type{<:AssemblyType} = ON_FACES, + xgrid = FE.dofgrid; + bonus_quadorder = 0, + operator = NormalFlux, + nfluxes = 0, + dofs = [], + mean = false, + kwargs... + ) where {Tv, Ti, FEType <: AbstractFiniteElement, APT} + itemvolumes = xgrid[GridComponentVolumes4AssemblyType(AT)] itemnodes = xgrid[GridComponentNodes4AssemblyType(AT)] itemregions = xgrid[GridComponentRegions4AssemblyType(AT)] @@ -493,7 +498,7 @@ function FunctionalInterpolator( end end if isempty(dofs) - dofs = 1 : nfluxes + dofs = 1:nfluxes end ncomponents::Int = get_ncomponents(FEType) order_FE = get_polynomialorder(FEType, EG) @@ -516,7 +521,7 @@ function FunctionalInterpolator( if item < 1 continue end - for m = 1:nfluxes + for m in 1:nfluxes target[itemdofs[dofs[m], item]] = 0 end QP.region = itemregions[item] @@ -529,7 +534,7 @@ function FunctionalInterpolator( if interior_offset > 0 update_basis!(FEB, item) end - + ## compute fluxes of function for qp in 1:nweights fill!(f_fluxes, 0) @@ -540,8 +545,8 @@ function FunctionalInterpolator( ## subtract flux of fixed dofs if interior_offset > 0 - for m = 1 : nfluxes, dof = 1 : interior_offset - f_fluxes[m] -= target[itemdofs[dof, item]] * FEB.cvals[m, dof, qp] + for m in 1:nfluxes, dof in 1:interior_offset + f_fluxes[m] -= target[itemdofs[dof, item]] * FEB.cvals[m, dof, qp] end end @@ -551,7 +556,7 @@ function FunctionalInterpolator( end ## set fluxes to dofs - for m = 1:nfluxes + for m in 1:nfluxes target[itemdofs[dofs[m], item]] += f_fluxes[m] * weight end end @@ -579,7 +584,6 @@ function FunctionalInterpolator( end - function slice(VTA::VariableTargetAdjacency, items = [], only_unique::Bool = true) subitems = zeros(Int, 0) if items == [] diff --git a/test/test_interpolators.jl b/test/test_interpolators.jl index 366fe5ec..d3251f9c 100644 --- a/test/test_interpolators.jl +++ b/test/test_interpolators.jl @@ -8,7 +8,7 @@ function run_interpolator_tests() H1P3{1, 1} => 3, H1Pk{1, 1, 3} => 3, H1Pk{1, 1, 4} => 4, - H1Pk{1, 1, 5} => 5 + H1Pk{1, 1, 5} => 5, ] TestCatalog2D = [ @@ -37,7 +37,7 @@ function run_interpolator_tests() H1P3{2, 2} => 3, H1Pk{2, 2, 3} => 3, H1Pk{2, 2, 4} => 4, - H1Pk{2, 2, 5} => 5 + H1Pk{2, 2, 5} => 5, ] TestCatalog3D = [ @@ -53,7 +53,7 @@ function run_interpolator_tests() H1P1TEB{3} => 1, H1BR{3} => 1, H1P2{3, 3} => 2, - H1P3{3, 3} => 3 + H1P3{3, 3} => 3, ] ## function that computes errors at enough quadrature points for polynomial of degree order diff --git a/test/test_segmentintegrator.jl b/test/test_segmentintegrator.jl index 4d70316e..08a9b2ad 100644 --- a/test/test_segmentintegrator.jl +++ b/test/test_segmentintegrator.jl @@ -37,7 +37,7 @@ function test_segmentintegrator_nokernel() error2 = sqrt((result[1] - 1 // 12)^2 + result[2]^2) println("error2 without kernel = $error2 (result = $result)") - return max(error1, error2) < 1e-15 + return max(error1, error2) < 1.0e-15 end function test_segmentintegrator_withkernel()