diff --git a/Project.toml b/Project.toml index 886c69b..9789446 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FiniteElementContainers" uuid = "d08262e4-672f-4e7f-a976-f2cea5767631" authors = ["Craig M. Hamel and contributors"] -version = "0.5.1" +version = "0.5.2" [deps] Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" @@ -28,7 +28,7 @@ Adapt = "3" Aqua = "0.8" Atomix = "0.1.0" DocStringExtensions = "0.9" -Exodus = "0.12" +Exodus = "0.13" IterativeSolvers = "0.9" JET = "0.8" KernelAbstractions = "0.9" diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index c39837b..013dd59 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -122,3 +122,54 @@ end include("DynamicAssemblers.jl") include("MatrixFreeStaticAssembler.jl") include("StaticAssembler.jl") + +""" +$(TYPEDSIGNATURES) +method that assumes first dof +TODO move sorting of nodes up stream +TODO remove other scratch unknowns and unknown_dofs arrays +""" +function update_unknown_dofs!( + assembler::Union{DynamicAssembler, StaticAssembler}, + dof, + fspaces, + nodes_in::V +) where V <: AbstractVector{<:Integer} + + # make this an assumption of the method + nodes = sort(nodes_in) + + n_total_dofs = num_dofs_per_node(dof) * num_nodes(dof) - length(nodes) + + # TODO change to a good sizehint! + resize!(assembler.Is, 0) + resize!(assembler.Js, 0) + resize!(assembler.unknown_dofs, 0) + + n = 1 + for fspace in fspaces + for e in 1:num_elements(fspace) + conn = dof_connectivity(fspace, e) + for temp in Iterators.product(conn, conn) + if insorted(temp[1], nodes) || insorted(temp[2], nodes) + # really do nothing here + else + push!(assembler.Is, temp[1] - count(x -> x < temp[1], nodes)) + push!(assembler.Js, temp[2] - count(x -> x < temp[2], nodes)) + push!(assembler.unknown_dofs, n) + end + n += 1 + end + end + end + + # resize cache arrays + resize!(assembler.klasttouch, n_total_dofs) + resize!(assembler.csrrowptr, n_total_dofs + 1) + resize!(assembler.csrcolval, length(assembler.Is)) + resize!(assembler.csrnzval, length(assembler.Is)) + + # resize!(assembler.stiffnesses, length(assembler.Is)) + + return nothing +end \ No newline at end of file diff --git a/src/assemblers/DynamicAssemblers.jl b/src/assemblers/DynamicAssemblers.jl index 0cb84b6..21b9e8d 100644 --- a/src/assemblers/DynamicAssemblers.jl +++ b/src/assemblers/DynamicAssemblers.jl @@ -12,6 +12,7 @@ struct DynamicAssembler{ U <: AbstractArray{Itype, 1}, Sizes <: AbstractArray{Itype, 1}, Offsets <: AbstractArray{Itype, 1}, + R <: NodalField, K <: AbstractArray{Rtype, 1}, # should always be a vector type thing M <: AbstractArray{Rtype, 1}, # cache types @@ -24,6 +25,7 @@ struct DynamicAssembler{ unknown_dofs::U block_sizes::Sizes block_offsets::Offsets + residuals::R stiffnesses::K masses::M # cache arrays @@ -78,6 +80,7 @@ function DynamicAssembler(dof::DofManager, fspaces::Fs) where Fs end end + residuals = create_fields(dof) #|> vec stiffnesses = zeros(Float64, size(Is)) masses = zeros(Float64, size(Is)) @@ -94,13 +97,15 @@ function DynamicAssembler(dof::DofManager, fspaces::Fs) where Fs return DynamicAssembler{ Float64, Int64, typeof(Is), typeof(Js), typeof(unknown_dofs), typeof(block_sizes), typeof(block_offsets), - typeof(stiffnesses), typeof(masses), + typeof(residuals), typeof(stiffnesses), typeof(masses), # cache arrays typeof(klasttouch), typeof(csrrowptr), typeof(csrcolval), typeof(csrnzval), # additional cache arrays typeof(csccolptr), typeof(cscrowval), typeof(cscnzval) }( - Is, Js, unknown_dofs, block_sizes, block_offsets, stiffnesses, masses, + Is, Js, unknown_dofs, block_sizes, block_offsets, + # fields + residuals, stiffnesses, masses, # cache arrays klasttouch, csrrowptr, csrcolval, csrnzval, # additional cache arrays @@ -163,7 +168,7 @@ $(TYPEDSIGNATURES) Simple method for assembling in serial """ function assemble!( - R, + # R, assembler::DynamicAssembler, dof::DofManager, fspace::FunctionSpace, @@ -175,6 +180,8 @@ function assemble!( N = num_nodes_per_element(fspace) NxNDof = N * NDof + R = assembler.residuals + for e in 1:num_elements(fspace) U_el = element_level_fields(fspace, U, e) R_el = zeros(SVector{NxNDof, Float64}) @@ -203,7 +210,7 @@ end $(TYPEDSIGNATURES) """ function assemble!( - R, + # R, assembler::DynamicAssembler, dof::DofManager, fspaces, X, U, @@ -211,11 +218,12 @@ function assemble!( ) # reset in some way + R = assembler.residuals R .= zero(eltype(R)) assembler.stiffnesses .= zero(eltype(assembler.stiffnesses)) for (block_id, fspace) in enumerate(fspaces) - assemble!(R, assembler, dof, fspace, X, U, block_id, residual_func, tangent_func, mass_func) + assemble!(assembler, dof, fspace, X, U, block_id, residual_func, tangent_func, mass_func) end end diff --git a/src/assemblers/StaticAssembler.jl b/src/assemblers/StaticAssembler.jl index 697b470..6ffe8f2 100644 --- a/src/assemblers/StaticAssembler.jl +++ b/src/assemblers/StaticAssembler.jl @@ -11,6 +11,7 @@ struct StaticAssembler{ U <: AbstractArray{Itype, 1}, Sizes <: AbstractArray{Itype, 1}, Offsets <: AbstractArray{Itype, 1}, + R <: NodalField, K <: AbstractArray{Rtype, 1}, # should always be a vector type thing # cache types C1, C2, C3, C4, @@ -22,6 +23,7 @@ struct StaticAssembler{ unknown_dofs::U block_sizes::Sizes block_offsets::Offsets + residuals::R stiffnesses::K # cache arrays klasttouch::C1 @@ -83,7 +85,7 @@ function StaticAssembler(dof::DofManager, fspaces::Fs) where Fs end end - # residuals = create_fields(dof) + residuals = create_fields(dof) #|> vec stiffnesses = zeros(Float64, size(Is)) # create caches @@ -99,13 +101,15 @@ function StaticAssembler(dof::DofManager, fspaces::Fs) where Fs return StaticAssembler{ Float64, Int64, typeof(Is), typeof(Js), typeof(unknown_dofs), typeof(block_sizes), typeof(block_offsets), - typeof(stiffnesses), + typeof(residuals), typeof(stiffnesses), # cache arrays typeof(klasttouch), typeof(csrrowptr), typeof(csrcolval), typeof(csrnzval), # additional cache arrays typeof(csccolptr), typeof(cscrowval), typeof(cscnzval) }( - Is, Js, unknown_dofs, block_sizes, block_offsets, stiffnesses, + Is, Js, unknown_dofs, block_sizes, block_offsets, + # fields + residuals, stiffnesses, # cache arrays klasttouch, csrrowptr, csrcolval, csrnzval, # additional cache arrays @@ -193,6 +197,8 @@ function assemble!( N = num_nodes_per_element(fspace) NxNDof = N * NDof + R = assembler.residuals + for e in 1:num_elements(fspace) U_el = element_level_fields(fspace, U, e) R_el = zeros(SVector{NxNDof, Float64}) @@ -221,7 +227,7 @@ Top level method using methods """ # TODO figure this one out function assemble!( - R, + # R, assembler::StaticAssembler, dof::DofManager, fspaces, X, U, @@ -230,6 +236,7 @@ function assemble!( # reset in some way # assembler.residuals .= 0. + R = assembler.residuals R .= zero(eltype(R)) assembler.stiffnesses .= zero(eltype(assembler.stiffnesses)) @@ -239,53 +246,53 @@ function assemble!( end -""" -$(TYPEDSIGNATURES) -method that assumes first dof -TODO move sorting of nodes up stream -TODO remove other scratch unknowns and unknown_dofs arrays -""" -function update_unknown_dofs!( - assembler::StaticAssembler, - dof, - fspaces, - nodes_in::V -) where V <: AbstractVector{<:Integer} - - # make this an assumption of the method - nodes = sort(nodes_in) - - n_total_dofs = num_dofs_per_node(dof) * num_nodes(dof) - length(nodes) - - # TODO change to a good sizehint! - resize!(assembler.Is, 0) - resize!(assembler.Js, 0) - resize!(assembler.unknown_dofs, 0) - - n = 1 - for fspace in fspaces - for e in 1:num_elements(fspace) - conn = dof_connectivity(fspace, e) - for temp in Iterators.product(conn, conn) - if insorted(temp[1], nodes) || insorted(temp[2], nodes) - # really do nothing here - else - push!(assembler.Is, temp[1] - count(x -> x < temp[1], nodes)) - push!(assembler.Js, temp[2] - count(x -> x < temp[2], nodes)) - push!(assembler.unknown_dofs, n) - end - n += 1 - end - end - end - - # resize cache arrays - resize!(assembler.klasttouch, n_total_dofs) - resize!(assembler.csrrowptr, n_total_dofs + 1) - resize!(assembler.csrcolval, length(assembler.Is)) - resize!(assembler.csrnzval, length(assembler.Is)) - - # resize!(assembler.stiffnesses, length(assembler.Is)) - - return nothing -end +# """ +# $(TYPEDSIGNATURES) +# method that assumes first dof +# TODO move sorting of nodes up stream +# TODO remove other scratch unknowns and unknown_dofs arrays +# """ +# function update_unknown_dofs!( +# assembler::StaticAssembler, +# dof, +# fspaces, +# nodes_in::V +# ) where V <: AbstractVector{<:Integer} + +# # make this an assumption of the method +# nodes = sort(nodes_in) + +# n_total_dofs = num_dofs_per_node(dof) * num_nodes(dof) - length(nodes) + +# # TODO change to a good sizehint! +# resize!(assembler.Is, 0) +# resize!(assembler.Js, 0) +# resize!(assembler.unknown_dofs, 0) + +# n = 1 +# for fspace in fspaces +# for e in 1:num_elements(fspace) +# conn = dof_connectivity(fspace, e) +# for temp in Iterators.product(conn, conn) +# if insorted(temp[1], nodes) || insorted(temp[2], nodes) +# # really do nothing here +# else +# push!(assembler.Is, temp[1] - count(x -> x < temp[1], nodes)) +# push!(assembler.Js, temp[2] - count(x -> x < temp[2], nodes)) +# push!(assembler.unknown_dofs, n) +# end +# n += 1 +# end +# end +# end + +# # resize cache arrays +# resize!(assembler.klasttouch, n_total_dofs) +# resize!(assembler.csrrowptr, n_total_dofs + 1) +# resize!(assembler.csrcolval, length(assembler.Is)) +# resize!(assembler.csrnzval, length(assembler.Is)) + +# # resize!(assembler.stiffnesses, length(assembler.Is)) + +# return nothing +# end diff --git a/test/poisson/TestPoisson.jl b/test/poisson/TestPoisson.jl index 98267b9..5636b10 100644 --- a/test/poisson/TestPoisson.jl +++ b/test/poisson/TestPoisson.jl @@ -58,9 +58,10 @@ for type in types function solve(asm, dof, fspaces, X, U, Uu) for n in 1:10 update_fields!(U, dof, Uu) - assemble!(R, asm, dof, fspaces, X, U, residual, tangent) + # assemble!(R, asm, dof, fspaces, X, U, residual, tangent) + assemble!(asm, dof, fspaces, X, U, residual, tangent) # R = asm.residuals[dof.unknown_dofs] - R_view = @views R[dof.unknown_dofs] + R_view = @views asm.residuals[dof.unknown_dofs] K = sparse(asm) cg!(ΔUu, -K, R_view) @show norm(ΔUu) norm(R_view)