diff --git a/Project.toml b/Project.toml index 99084eb..274fe61 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableFEMBase" uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c" authors = ["Christian Merdon "] -version = "0.3.1" +version = "0.3.2" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" diff --git a/examples/Example205_LowLevelSpaceTimePoisson.jl b/examples/Example205_LowLevelSpaceTimePoisson.jl new file mode 100644 index 0000000..7fb0841 --- /dev/null +++ b/examples/Example205_LowLevelSpaceTimePoisson.jl @@ -0,0 +1,317 @@ +#= + +# 205 : Space-Time FEM for Poisson Problem +([source code](SOURCE_URL)) + +This example computes the solution ``u`` of the +two-dimensional heat equation +```math +\begin{aligned} +u_t - \Delta u & = f \quad \text{in } \Omega +\end{aligned} +``` +with a (possibly space- and time-depepdent) +right-hand side ``f`` and homogeneous Dirichlet boundary +and initial conditions +on the unit square domain ``\Omega`` on a given grid +with space-time finite element methods based on +tensorized ansatz functions. + +![](example205.svg) + +=# + +module Example205_LowLevelSpaceTimePoisson + +using ExtendableFEMBase +using ExtendableGrids +using ExtendableSparse +using GridVisualize +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) + +function main(; dt = 0.01, Tfinal = 1, level = 5, order = 1, produce_movie = true, Plotter = nothing) + + ## 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)) + grid_time = simplexgrid(T) + + ## space grid + 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 + 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 + +end + + +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 + bfacedofs::Adjacency{Int32} = FES_space[ExtendableFEMBase.BFaceDofs] + nbfaces::Int = num_sources(bfacedofs) + dof::Int = 0 + for bface ∈ 1:nbfaces + for dof_t = 1 : ndofs_time + for j ∈ 1:num_targets(bfacedofs, 1) + dof = (dof_t-1)*ndofs_space + bfacedofs[j, bface] + A[dof, dof] = 1e60 + b[dof] = 0 + end + end + end + + ## fix initial value by zero + for j=1:ndofs_space + A[j,j] = 1e60 + b[dof] = 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) + + ## compute linear residual + @show norm(A*sol - b) + + 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 + 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) + Mloc = zeros(Float64, ndofs4cell_time, ndofs4cell_time) + 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) + + ## 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 + 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 + dof_j = celldofs_space[j, cell] + (dofTj - 1) * ndofs_space + for k ∈ 1:ndofs4cell_space + dof_k = celldofs_space[k, cell] + (dofTk - 1) * ndofs_space + if abs(Aloc[j, k]) > 1e-15 + ## write into sparse matrix, only lines with allocations + rawupdateindex!(A, +, Aloc[j, k]*factor, dof_j, dof_k) + 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 + 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 + dofXj = celldofs_space[jX, cell] + dofXk = celldofs_space[kX, cell] + for qpX ∈ 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 + dof_j = dofXj + (celldofs_time[j, time_cell] - 1) * ndofs_space + for k ∈ 1:ndofs4cell_time + dof_k = dofXk + (celldofs_time[k, time_cell] - 1) * ndofs_space + if abs(Mloc[j, k]) > 1e-15 + ## write into sparse matrix, only lines with allocations + rawupdateindex!(A, +, Mloc[j, k]*factor, dof_j, dof_k) + 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 +end + +function generateplots(dir = pwd(); Plotter = nothing, kwargs...) + ~, plt = main(; Plotter = Plotter, kwargs...) + scene = GridVisualize.reveal(plt) + GridVisualize.save(joinpath(dir, "example205.svg"), scene; Plotter = Plotter) +end + +end #module \ No newline at end of file diff --git a/src/feevaluator.jl b/src/feevaluator.jl index 013725d..7f1c371 100644 --- a/src/feevaluator.jl +++ b/src/feevaluator.jl @@ -179,7 +179,7 @@ _prepare_additional_coefficients(::Type{<:NormalFlux}, xgrid, AT, edim) = AT == _prepare_additional_coefficients(::Type{<:TangentialGradient}, xgrid, AT, edim) = xgrid[FaceNormals] function _prepare_additional_coefficients(::Type{<:TangentFlux}, xgrid, AT, edim) if edim == 2 - return _prepare_additional_coefficients(NormalFlux, AT, edim) + return _prepare_additional_coefficients(NormalFlux, xgrid, AT, edim) else return AT == ON_BEDGES ? view(xgrid[EdgeTangents], :, xgrid[BEdgeEdges]) : xgrid[EdgeTangents] end diff --git a/src/feevaluator_hcurl.jl b/src/feevaluator_hcurl.jl index 7a0b9fc..7383ed4 100644 --- a/src/feevaluator_hcurl.jl +++ b/src/feevaluator_hcurl.jl @@ -22,6 +22,7 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Tang subset = _update_subset!(FEBE) refbasisvals = FEBE.refbasisvals xItemVolumes = FEBE.L2G.ItemVolumes + item = FEBE.citem[] cvals = FEBE.cvals for i ∈ 1:size(cvals, 3), dof_i ∈ 1:size(cvals, 2), k ∈ 1:size(cvals, 1) cvals[k, dof_i, i] = refbasisvals[i][subset[dof_i], k] / xItemVolumes[item] diff --git a/src/finiteelements.jl b/src/finiteelements.jl index 67b92c9..08dc019 100644 --- a/src/finiteelements.jl +++ b/src/finiteelements.jl @@ -95,11 +95,14 @@ function FESpace{FEType, AT}( xgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 0) xgrid[BFaceVolumes] = zeros(Tv, 0) elseif AT == ON_BFACES - xgrid = subgrid(xgrid, unique(xgrid[BFaceRegions]); boundary = true, project = false) - xgrid[BFaceNodes] = zeros(Ti, 2, 0) - xgrid[BFaceRegions] = zeros(Ti, 0) - xgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 0) - xgrid[BFaceVolumes] = zeros(Tv, 0) + if isnothing(regions) + regions = unique(xgrid[BFaceRegions]) + end + dofgrid = subgrid(xgrid, regions; boundary = true, project = false) + dofgrid[BFaceNodes] = zeros(Ti, 2, 0) + dofgrid[BFaceRegions] = zeros(Ti, 0) + dofgrid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Int}(Edge1D, 0) + dofgrid[BFaceVolumes] = zeros(Tv, 0) elseif AT == ON_EDGES xgrid = get_edgegrid(xgrid) xgrid[BFaceNodes] = zeros(Ti, 2, 0) @@ -112,10 +115,12 @@ function FESpace{FEType, AT}( regions = unique(xgrid[CellRegions]) end - if regions != unique(xgrid[CellRegions]) - dofgrid = subgrid(xgrid, regions) - else - dofgrid = xgrid + if AT !== ON_BFACES + if regions != unique(xgrid[CellRegions]) + dofgrid = subgrid(xgrid, regions) + else + dofgrid = xgrid + end end # first generate some empty FESpace diff --git a/src/functionoperators.jl b/src/functionoperators.jl index ecd3db8..9b33558 100644 --- a/src/functionoperators.jl +++ b/src/functionoperators.jl @@ -168,7 +168,7 @@ Length4Operator(::Type{CurlScalar}, xdim::Int, ncomponents::Int) = ((xdim == 2) Length4Operator(::Type{Curl2D}, xdim::Int, ncomponents::Int) = 1 Length4Operator(::Type{Curl3D}, xdim::Int, ncomponents::Int) = 3 Length4Operator(::Type{<:Gradient}, xdim::Int, ncomponents::Int) = xdim * ncomponents -Length4Operator(::Type{TangentialGradient}, xdim::Int, ncomponents::Int) = 1 +Length4Operator(::Type{TangentialGradient}, xdim::Int, ncomponents::Int) = ncomponents Length4Operator(::Type{<:SymmetricGradient}, xdim::Int, ncomponents::Int) = ((xdim == 2) ? 3 : 6) * Int(ceil(ncomponents / xdim)) Length4Operator(::Type{<:Hessian}, xdim::Int, ncomponents::Int) = xdim * xdim * ncomponents Length4Operator(::Type{<:SymmetricHessian}, xdim::Int, ncomponents::Int) = ((xdim == 2) ? 3 : 6) * ncomponents diff --git a/src/interpolations.jl b/src/interpolations.jl index bad46b9..56fe864 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -1042,22 +1042,19 @@ function displace_mesh!(xgrid::ExtendableGrid, source::FEVectorBlock; magnify = xgrid[Coordinates] .+= magnify * nodevals # remove all keys from grid components that might have changed and need a reinstantiation - delete!(xgrid.components, CellVolumes) - delete!(xgrid.components, FaceVolumes) - delete!(xgrid.components, EdgeVolumes) - delete!(xgrid.components, FaceNormals) - delete!(xgrid.components, EdgeTangents) - delete!(xgrid.components, BFaceVolumes) - delete!(xgrid.components, BEdgeVolumes) + ExtendableGrids.update!(xgrid, CellVolumes) + ExtendableGrids.update!(xgrid, FaceVolumes) + ExtendableGrids.update!(xgrid, EdgeVolumes) + ExtendableGrids.update!(xgrid, FaceNormals) + ExtendableGrids.update!(xgrid, EdgeTangents) + ExtendableGrids.update!(xgrid, BFaceVolumes) + ExtendableGrids.update!(xgrid, BEdgeVolumes) end -function displace_mesh(xgrid::ExtendableGrid, source::FEVectorBlock; magnify = 1) +function displace_mesh(xgrid::ExtendableGrid, source::FEVectorBlock; kwargs...) xgrid_displaced = deepcopy(xgrid) - nnodes = size(xgrid[Coordinates], 2) - nodevals = zeros(eltype(xgrid[Coordinates]), get_ncomponents(Base.eltype(source.FES)), nnodes) - nodevalues!(nodevals, source, Identity) - xgrid_displaced[Coordinates] .+= magnify * nodevals + displace_mesh!(xgrid_displaced, source; kwargs...) return xgrid_displaced end