Merge pull request #34 from Cthonios/dev

Dev
Cthonios · Jan 18, 2024 · 05dca9e · 05dca9e
2 parents 0196537 + 5ae3233
commit 05dca9e
Show file tree

Hide file tree

Showing 20 changed files with 910 additions and 734 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,11 +1,10 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
-version = "0.3.1"
+version = "0.4.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ReferenceFiniteElements = "6dc62d09-f8eb-43fd-9672-074e490a997f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -26,10 +25,8 @@ FiniteElementContainersKernelAbstractionsExt = "KernelAbstractions"
 [compat]
 Adapt = "3"
 Aqua = "0.8"
-CUDA = "5"
 DocStringExtensions = "0.9"
 Exodus = "0.11"
-FillArrays = "1"
 IterativeSolvers = "0.9"
 JET = "0.8"
 KernelAbstractions = "0.9"
@@ -46,12 +43,12 @@ julia = "1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 
 [targets]
-test = ["Aqua", "CUDA", "Exodus", "IterativeSolvers", "JET", "Parameters", "Test", "TestSetExtensions"]
+test = ["Aqua", "Exodus", "IterativeSolvers", "JET", "Parameters", "Test", "TestSetExtensions"]
diff --git a/ext/FiniteElementContainersAdaptExt.jl b/ext/FiniteElementContainersAdaptExt.jl
@@ -39,22 +39,9 @@ end
 # no need for adapt_structure for connectivity fields since these are aliases anyway
 
 # DofManagers
-function Adapt.adapt_structure(to, dof::FiniteElementContainers.SimpleDofManager)
-  ND, NN = num_dofs_per_node(dof), num_nodes(dof)
-  is_unknown      = Adapt.adapt_structure(to, dof.is_unknown)
-  unknown_indices = Adapt.adapt_structure(to, dof.unknown_indices)
-  return FiniteElementContainers.SimpleDofManager{Bool, 2, ND, NN, typeof(is_unknown), typeof(unknown_indices)}(
-    is_unknown, unknown_indices
-  )
-end
-
-function Adapt.adapt_structure(to, dof::FiniteElementContainers.VectorizedDofManager)
-  ND, NN = num_dofs_per_node(dof), num_nodes(dof)
-  is_unknown      = Adapt.adapt_structure(to, dof.is_unknown)
-  unknown_indices = Adapt.adapt_structure(to, dof.unknown_indices)
-  return FiniteElementContainers.VectorizedDofManager{Bool, 2, ND, NN, typeof(is_unknown), typeof(unknown_indices)}(
-    is_unknown, unknown_indices
-  )
+function Adapt.adapt_structure(to, dof::FiniteElementContainers.DofManager{T, ND, NN, A, V}) where {T, ND, NN, A, V}
+  unknown_dofs = Adapt.adapt_structure(to, dof.unknown_dofs)
+  return FiniteElementContainers.DofManager{T, ND, NN, A, typeof(unknown_dofs)}(unknown_dofs)
 end
 
 # Function spaces
@@ -69,15 +56,57 @@ function Adapt.adapt_structure(to, fspace::FiniteElementContainers.NonAllocatedF
 end
 
 # Assemblers
-# """
-# Need to use SparseArrays.allowscalar(false)
-# """
-# function Adapt.adapt_structure(to, assembler::FiniteElementContainers.StaticAssembler)
-#   I = FiniteElementContainers.int_type(assembler)
-#   F = FiniteElementContainers.float_type(assembler)
-#   R = Adapt.adapt_structure(to, assembler.R)
-#   K = Adapt.adapt_structure(to, assembler.K)
-#   return FiniteElementContainers.StaticAssembler{I, F, typeof(R), typeof(K)}(R, K)
-# end
+function Adapt.adapt_structure(to, asm::FiniteElementContainers.DynamicAssembler)
+  I = FiniteElementContainers.int_type(assembler)
+  F = FiniteElementContainers.float_type(assembler)
+
+  Is = Adapt.adapt_structure(to, asm.Is)
+  Js = Adapt.adapt_structure(to, asm.Js)
+  unknown_dofs = Adapt.adapt_structure(to, asm.unknown_dofs)
+  block_sizes = Adapt.adapt_structure(to, asm.block_sizes)
+  block_offsets = Adapt.adapt_structure(to, asm.block_offsets)
+  residuals = Adapt.adapt_structure(to, asm.residuals)
+  stiffnesses = Adapt.adapt_structure(to, asm.stiffnesses)
+  masses = Adapt.adapt_structure(to, asm.masses)
+
+  return FiniteElementContainers.StaticAssembler{
+    F, I,
+    typeof(Is), typeof(Js), 
+    typeof(unknown_dofs), typeof(block_sizes), typeof(block_offsets),
+    typeof(residuals), tyepof(stiffnesses), typeof(masses)
+  }(Is, Js, unknown_dofs, block_sizes, block_offsets, residuals, stiffnesses, masses)
+end
+
+function Adapt.adapt_structure(to, asm::FiniteElementContainers.StaticAssembler)
+  I = FiniteElementContainers.int_type(assembler)
+  F = FiniteElementContainers.float_type(assembler)
+
+  Is = Adapt.adapt_structure(to, asm.Is)
+  Js = Adapt.adapt_structure(to, asm.Js)
+  unknown_dofs = Adapt.adapt_structure(to, asm.unknown_dofs)
+  block_sizes = Adapt.adapt_structure(to, asm.block_sizes)
+  block_offsets = Adapt.adapt_structure(to, asm.block_offsets)
+  residuals = Adapt.adapt_structure(to, asm.residuals)
+  stiffnesses = Adapt.adapt_structure(to, asm.stiffnesses)
+
+  return FiniteElementContainers.StaticAssembler{
+    F, I,
+    typeof(Is), typeof(Js), 
+    typeof(unknown_dofs), typeof(block_sizes), typeof(block_offsets),
+    typeof(residuals), tyepof(stiffnesses)
+  }(Is, Js, unknown_dofs, block_sizes, block_offsets, residuals, stiffnesses)
+end 
+
+
+# # """
+# # Need to use SparseArrays.allowscalar(false)
+# # """
+# # function Adapt.adapt_structure(to, assembler::FiniteElementContainers.StaticAssembler)
+# #   I = FiniteElementContainers.int_type(assembler)
+# #   F = FiniteElementContainers.float_type(assembler)
+# #   R = Adapt.adapt_structure(to, assembler.R)
+# #   K = Adapt.adapt_structure(to, assembler.K)
+# #   return FiniteElementContainers.StaticAssembler{I, F, typeof(R), typeof(K)}(R, K)
+# # end
 
 end # module
diff --git a/ext/FiniteElementContainersExodusExt.jl b/ext/FiniteElementContainersExodusExt.jl
@@ -3,18 +3,22 @@ module FiniteElementContainersExodusExt
 using Exodus
 using FiniteElementContainers
 
-function FiniteElementContainers.FileMesh(::Type{ExodusDatabase}, file_name::String)
-  exo = ExodusDatabase(file_name, "r")
+function FiniteElementContainers.FileMesh(type::Type{<:ExodusDatabase}, file_name::String)
+  exo = type(file_name, "r")
   return FileMesh{typeof(exo)}(file_name, exo)
 end
 
-FiniteElementContainers.num_dimensions(
+function FiniteElementContainers.num_dimensions(
   mesh::FileMesh{<:ExodusDatabase}
-)  = mesh.mesh_obj.init.num_dim
+)::Int32
+  return mesh.mesh_obj.init.num_dim
+end
 
-FiniteElementContainers.num_nodes(
+function FiniteElementContainers.num_nodes(
   mesh::FileMesh{<:ExodusDatabase}
-) = mesh.mesh_obj.init.num_nodes
+)::Int32
+  return mesh.mesh_obj.init.num_nodes
+end
 
 function FiniteElementContainers.element_block_ids(mesh::FileMesh{<:ExodusDatabase})
   return Exodus.read_ids(mesh.mesh_obj, Block)
@@ -28,7 +32,7 @@ function FiniteElementContainers.sideset_ids(mesh::FileMesh{<:ExodusDatabase})
   return Exodus.read_ids(mesh.mesh_obj, SideSet)
 end
 
-function FiniteElementContainers.coordinates(mesh::FileMesh{<:ExodusDatabase}) 
+function FiniteElementContainers.coordinates(mesh::FileMesh{ExodusDatabase{M, I, B, F}})::Matrix{F} where {M, I, B, F} 
   coords = Exodus.read_coordinates(mesh.mesh_obj)
   return coords
 end
@@ -70,4 +74,4 @@ function FiniteElementContainers.sideset(
   return sset.elements, sset.sides
 end
 
-end # module
+end # module
diff --git a/ext/FiniteElementContainersKernelAbstractionsExt.jl b/ext/FiniteElementContainersKernelAbstractionsExt.jl
@@ -5,72 +5,106 @@ using KernelAbstractions
 
 # DofManagers
 function FiniteElementContainers.create_fields(
-  backend::Backend, 
-  ::DofManager{T, N, ND, NN, B}, 
-  float_type::Type{<:Number} = Float64
-) where {T, N, ND, NN, B <: AbstractArray{Bool, 1}}
-  vals = KernelAbstractions.zeros(backend, float_type, ND * NN)
+  ::DofManager{T, ND, NN, A, V},
+  backend::Backend
+) where {T, ND, NN, A <: AbstractVector, V}
+  vals = KernelAbstractions.zeros(backend, eltype(A), ND * NN)
   return NodalField{ND, NN, Vector}(vals)
 end
 
 function FiniteElementContainers.create_fields(
-  backend::Backend, 
-  ::DofManager{T, N, ND, NN, B}, 
-  float_type::Type{<:Number} = Float64
-) where {T, N, ND, NN, B <: AbstractArray{Bool, 2}}
-  vals = KernelAbstractions.zeros(backend, float_type, ND, NN)
+  ::DofManager{T, ND, NN, A, V},
+  backend::Backend
+) where {T, ND, NN, A <: AbstractMatrix, V}
+  vals = KernelAbstractions.zeros(backend, eltype(A), ND, NN)
   return NodalField{ND, NN, Matrix}(vals)
 end
 
 function FiniteElementContainers.create_unknowns(
-  backend::Backend,
-  dof::DofManager{T, N, ND, NN, Bools}, 
-  float_type::Type{<:Number} = Float64
-) where {T, N, ND, NN, Bools}
-  return KernelAbstractions.zeros(backend, float_type, length(dof.unknown_indices))
+  dof::DofManager{T, ND, NN, A, V},
+  backend::Backend
+) where {T, ND, NN, A, V}
+  return KernelAbstractions.zeros(backend, eltype(A), length(dof.unknown_dofs))
 end
 
 @kernel function update_fields_kernel!(U::NodalField, dof::DofManager, Uu::V) where V <: AbstractVector
   I = @index(Global)
-  U[dof.unknown_indices[I]] = Uu[I]
+  U[dof.unknown_dofs[I]] = Uu[I]
 end
 
 # TODO add wg size as parameter
-function FiniteElementContainers.update_fields!(backend::Backend, U::NodalField, dof::DofManager, Uu::V) where V <: AbstractVector
-  @assert length(Uu) == sum(dof.is_unknown)
+function FiniteElementContainers.update_fields!(U::NodalField, dof::DofManager, Uu::V, backend::Backend) where V <: AbstractVector
+  @assert length(Uu) == length(dof.unknown_dofs)#sum(dof.is_unknown)
   kernel = update_fields_kernel!(backend)
   kernel(U, dof, Uu, ndrange=length(Uu))
   synchronize(backend)
 end
 
-@kernel function update_unknown_ids_nset_kernel!(dof_manager, nset, dof)
+@kernel function update_unknown_dofs_kernel!(dof, dofs)
   I = @index(Global)
-  dof_manager[dof, nset[I]] = 0
+  dof.unknown_dofs[I] = dofs[I]
+  nothing
 end
 
-@kernel function update_unknown_ids_kernel!(dof, ids)
-  I = @index(Global)
-  dof.unknown_indices[I] = ids[dof.unknown_indices[I]]
+function FiniteElementContainers.update_unknown_dofs!(
+  dof::DofManager{T, ND, NN, A, V},
+  backend::Backend 
+) where {T, ND, NN, A, V}
+  resize!(dof.unknown_dofs, ND * NN)
+  dofs = FiniteElementContainers.dof_ids(dof)
+  kernel = update_unknown_dofs_kernel!(backend)
+  kernel(dof, dofs, ndrange=length(dofs))
+  nothing
 end
 
-function FiniteElementContainers.update_unknown_ids!(
-  backend::Backend,
-  dof::DofManager,
-  nsets::Nsets,
-  dofs::Dofs
-) where {Nsets, Dofs}
+# function FiniteElementContainers.update_unknown_dofs!(
+#   dof::DofManager{T, ND, NN, A, V1},
+#   dofs::V2,
+#   backend::Backend 
+# ) where {T, ND, NN, A, V1, V2 <: AbstractVector{<:Integer}}
+#   resize!(dof.unknown_dofs, ND * NN)
+#   dofs_all = FiniteElementContainers.dof_ids(dof)
+#   kernel = update_unknown_dofs_kernel!(backend)
+#   kernel(dof, dofs_all, ndrange=length(dofs_all))
+#   # setdiff!(dofs_all, dofs)
+#   # resize!(dof.unknown_indices, length(dofs_all))
+#   # dof.unknown_indices .= dofs_all
+#   # deleteat!(dof.unknown_dofs, dofs)
+#   setdiff!(dof.unknown_dofs, dofs)
+#   nothing
+# end
 
-  kernel1 = update_unknown_ids_nset_kernel!(backend)
-  kernel2 = update_unknown_ids_kernel!(backend)
-  @assert length(nsets) == length(dofs)
-  for (n, nset) in enumerate(nsets)
-    @assert dofs[n] > 0 && dofs[n] <= num_dofs_per_node(dof)
-    kernel1(dof, nset, dofs[n], ndrange=length(nset))
-  end
-  resize!(dof.unknown_indices, sum(dof.is_unknown))
-  ids = FiniteElementContainers.dof_ids(dof)
-  kernel2(dof, ids, ndrange=length(dof.unknown_indices))
-  synchronize(backend)
-end
+# @kernel function update_unknown_dofs_non_default_kernel!(dof, dofs)
+#   I
+# end
+# @kernel function update_unknown_ids_nset_kernel!(dof_manager, nset, dof)
+#   I = @index(Global)
+#   dof_manager[dof, nset[I]] = 0
+# end
+
+# @kernel function update_unknown_ids_kernel!(dof, ids)
+#   I = @index(Global)
+#   dof.unknown_indices[I] = ids[dof.unknown_dofs[I]]
+# end
+
+# function FiniteElementContainers.update_unknown_dofs!(
+#   dof::DofManager,
+#   nsets::Nsets,
+#   dofs::Dofs,
+#   backend::Backend
+# ) where {Nsets, Dofs}
+
+#   kernel1 = update_unknown_ids_nset_kernel!(backend)
+#   kernel2 = update_unknown_ids_kernel!(backend)
+#   @assert length(nsets) == length(dofs)
+#   for (n, nset) in enumerate(nsets)
+#     @assert dofs[n] > 0 && dofs[n] <= num_dofs_per_node(dof)
+#     kernel1(dof, nset, dofs[n], ndrange=length(nset))
+#   end
+#   resize!(dof.unknown_indices, length(dof.unknown_dofs))
+#   ids = FiniteElementContainers.dof_ids(dof)
+#   kernel2(dof, ids, ndrange=length(dof.unknown_dofs))
+#   synchronize(backend)
+# end
 
 end # module