some more testing.

ext/FiniteElementContainersAdaptExt.jl

@@ -60,20 +60,23 @@ function Adapt.adapt_structure(to, bk::FiniteElementContainers.BCBookKeeping{V})
 end
 
 function Adapt.adapt_structure(to, bc::FiniteElementContainers.DirichletBCContainer)
-  bk = adapt(to, bc.bookkeeping)
+  dofs = adapt(to, bc.dofs)
+  nodes = adapt(to, bc.nodes)
   vals = adapt(to, bc.vals)
   vals_dot = adapt(to, bc.vals_dot)
   vals_dot_dot = adapt(to, bc.vals_dot_dot)
-  return FiniteElementContainers.DirichletBCContainer(bk, vals, vals_dot, vals_dot_dot)
+  return FiniteElementContainers.DirichletBCContainer(dofs, nodes, vals, vals_dot, vals_dot_dot)
 end
 
 function Adapt.adapt_structure(to, bc::FiniteElementContainers.NeumannBCContainer)
-  bk = adapt(to, bc.bookkeeping)
   el_conns = adapt(to, bc.element_conns)
+  elements = adapt(to, bc.elements)
+  side_nodes = adapt(to, bc.side_nodes)
+  sides = adapt(to, bc.sides)
   surf_conns = adapt(to, bc.surface_conns)
   ref_fe = adapt(to, bc.ref_fe)
   vals = adapt(to, bc.vals)
-  return FiniteElementContainers.NeumannBCContainer(bk, el_conns, surf_conns, ref_fe, vals)
+  return FiniteElementContainers.NeumannBCContainer(el_conns, elements, side_nodes, sides, surf_conns, ref_fe, vals)
 end
 
 # DofManagers

src/Parameters.jl

@@ -196,6 +196,12 @@ function create_parameters(
   return Parameters(mesh, assembler, physics, props, ics, dirichlet_bcs, neumann_bcs, times)
 end
 
+function initialize!(p::Parameters)
+  update_ic_values!(p.ics, p.ic_funcs, p.h1_coords)
+  update_field_ics!(p.h1_field, p.ics)
+  return nothing
+end
+
 """
 $(TYPEDSIGNATURES)
 This method is used to update the stored bc values.

src/assemblers/NeumannBC.jl

@@ -39,7 +39,7 @@ function _assemble_block_vector_neumann_bc!(
   for e in axes(conns, 2)
     x_el = _element_level_fields(X, ref_fe, conns, e)
     R_el = _element_scratch_vector(surface_element(ref_fe.element), U)
-    side = bc.bookkeeping.sides[e]
+    side = bc.sides[e]
     for q in 1:num_quadrature_points(surface_element(ref_fe.element))
       interps = MappedSurfaceInterpolants(ref_fe, x_el, q, side)
       Nvec = interps.N_reduced
@@ -54,7 +54,7 @@ function _assemble_block_vector_neumann_bc!(
       R_el = R_el + JxW * Nvec * f_val
     end
     block_id = 1 # doesn't matter for this method
-    el_id = bc.bookkeeping.elements[e]
+    el_id = bc.elements[e]
     @views _assemble_element!(field, R_el, bc.surface_conns[:, e], el_id, block_id)
   end
 end
@@ -83,7 +83,7 @@ KA.@kernel function _assemble_block_vector_neumann_bc_kernel!(
 
   x_el = _element_level_fields(X, ref_fe, conns, E)
   R_el = _element_scratch_vector(surface_element(ref_fe.element), U)
-  side = bc.bookkeeping.sides[E]
+  side = bc.sides[E]
 
   for q in 1:num_quadrature_points(surface_element(ref_fe.element))
     interps = MappedSurfaceInterpolants(ref_fe, x_el, q, side)

src/assemblers/SparseMatrixAssembler.jl

@@ -228,7 +228,7 @@ function update_dofs!(assembler::SparseMatrixAssembler, dirichlet_bcs)
   use_condensed = _is_condensed(assembler.dof)
 
   if length(dirichlet_bcs) > 0
-    dirichlet_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, dirichlet_bcs)
+    dirichlet_dofs = mapreduce(x -> x.dofs, vcat, dirichlet_bcs)
     dirichlet_dofs = unique(sort(dirichlet_dofs))
   else
     dirichlet_dofs = Vector{Int}(undef, 0)

src/bcs/BoundaryConditions.jl

@@ -28,6 +28,8 @@ function _unique_sort_perm(array::AbstractArray{T, 1}) where T <: Number
   return [id_map[x] for x in sorted_unique]
 end
 
+const SetName = Union{Nothing, Symbol}
+
 """
 $(TYPEDEF)
 $(TYPEDSIGNATURES)
@@ -55,36 +57,81 @@ end
 $(TYPEDSIGNATURES)
 """
 function BCBookKeeping(
-  mesh, dof::DofManager, var_name::Symbol, sset_name::Symbol
+  mesh, dof::DofManager, var_name::Symbol; #sset_name::Symbol
+  block_name::SetName = nothing,
+  nset_name::SetName = nothing,
+  sset_name::SetName = nothing
 )
+  # check to ensure at least one name is supplied
+  if block_name === nothing &&
+     nset_name === nothing &&
+     sset_name === nothing
+    @assert false "Need to specify either a block, nodeset or sideset."
+  end
+
+  # now check to make sure only one name is not nothing
+  if block_name !== nothing
+    @assert nset_name === nothing && sset_name === nothing
+  elseif nset_name !== nothing
+    @assert block_name === nothing && sset_name === nothing
+  elseif sset_name !== nothing  
+    @assert block_name === nothing && nset_name === nothing
+  end 
+
   # get dof index associated with this var
   dof_index = _dof_index_from_var_name(dof, var_name)
-
   fspace = function_space(dof)
 
   # get sset specific fields
-  elements = getproperty(mesh.sideset_elems, sset_name)
-  nodes = getproperty(mesh.sideset_nodes, sset_name)
-  sides = getproperty(mesh.sideset_sides, sset_name)
-  side_nodes = getproperty(mesh.sideset_side_nodes, sset_name)
-
-  blocks = Vector{Int64}(undef, 0)
-
-  # gather the blocks that are present in this sideset
-  # TODO this isn't quite right
-  for (n, val) in enumerate(values(mesh.element_id_maps))
-    # note these are the local elem id to the block, e.g. starting from 1.
-    indices_in_sset = indexin(val, elements)
-    filter!(x -> x !== nothing, indices_in_sset)
-
-    if length(indices_in_sset) > 0
-      append!(blocks, repeat([n], length(indices_in_sset)))
+  all_dofs = reshape(1:length(dof), size(dof))
+
+  if block_name !== nothing
+    # for this case it is likely a DirichletBC or InitialCondition
+    # so we really only need the nodes/dofs although this might
+    # not be the case for say Hdiv or Hcurl fields...
+    # TODO eventually set the blocks, could be useful maybe?
+    blocks = Vector{Int64}(undef, 0)
+    conns = getproperty(mesh.element_conns, block_name)
+    nodes = sort(unique(conns.data))
+    dofs = all_dofs[dof_index, nodes]
+    elements = getproperty(mesh.element_id_maps, block_name)
+    # below 2 don't make sense for other mesh entity types
+    sides = Vector{Int64}(undef, 0)
+    side_nodes = Matrix{Int64}(undef, 0, 0)
+  elseif nset_name !== nothing
+    # for this case we only setup the "nodes" and "dofs"
+    blocks = Vector{Int64}(undef, 0) # TODO we could eventually put the blocks present here
+    nodes = getproperty(mesh.nodeset_nodes, nset_name)
+    dofs = all_dofs[dof_index, nodes]
+    elements = Vector{Int64}(undef, 0) # TODO we could eventually put the elements present here
+    # below 2 don't make sense for other mesh entity types
+    sides = Vector{Int64}(undef, 0)
+    side_nodes = Matrix{Int64}(undef, 0, 0)
+  elseif sset_name !== nothing
+
+    elements = getproperty(mesh.sideset_elems, sset_name)
+    nodes = getproperty(mesh.sideset_nodes, sset_name)
+    sides = getproperty(mesh.sideset_sides, sset_name)
+    side_nodes = getproperty(mesh.sideset_side_nodes, sset_name)
+
+    blocks = Vector{Int64}(undef, 0)
+
+    # gather the blocks that are present in this sideset
+    # TODO this isn't quite right
+    for (n, val) in enumerate(values(mesh.element_id_maps))
+      # note these are the local elem id to the block, e.g. starting from 1.
+      indices_in_sset = indexin(val, elements)
+      filter!(x -> x !== nothing, indices_in_sset)
+
+      if length(indices_in_sset) > 0
+        append!(blocks, repeat([n], length(indices_in_sset)))
+      end
     end
-  end
 
-  # setup dofs local to this BC
-  all_dofs = reshape(1:length(dof), size(dof))
-  dofs = all_dofs[dof_index, nodes]
+    # setup dofs local to this BC
+    # all_dofs = reshape(1:length(dof), size(dof))
+    dofs = all_dofs[dof_index, nodes]
+  end
 
   return BCBookKeeping(
     blocks, dofs, elements, nodes, sides, side_nodes
@@ -143,14 +190,15 @@ $(TYPEDEF)
 $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
-abstract type AbstractBCContainer{
-  IT <: Integer,
-  VT <: Union{<:Number, <:SVector},
-  N,
-  IV <: AbstractArray{IT, 1},
-  IM <: AbstractArray{IT, 2},
-  VV <: AbstractArray{VT, N}
-} end
+# abstract type AbstractBCContainer{
+#   IT <: Integer,
+#   VT <: Union{<:Number, <:SVector},
+#   N,
+#   IV <: AbstractArray{IT, 1},
+#   IM <: AbstractArray{IT, 2},
+#   VV <: AbstractArray{VT, N}
+# } end
+abstract type AbstractBCContainer end
 
 KA.get_backend(x::AbstractBCContainer) = KA.get_backend(x.vals)
 

src/bcs/DirichletBCs.jl

@@ -1,7 +1,7 @@
 abstract type AbstractDirichletBC{F} <: AbstractBC{F} end
-abstract type AbstractDirichletBCContainer{
-  IT, VT, IV, IM, VV
-} <: AbstractBCContainer{IT, VT, 1, IV, IM, VV} end
+# abstract type AbstractDirichletBCContainer{
+#   IT, VT, IV, IM, VV
+# } <: AbstractBCContainer{IT, VT, 1, IV, IM, VV} end
 
 """
 $(TYPEDEF)
@@ -43,10 +43,10 @@ struct DirichletBCContainer{
   IT <: Integer,
   VT <: Number,
   IV <: AbstractArray{IT, 1},
-  IM <: AbstractArray{IT, 2},
   VV <: AbstractArray{VT, 1}
-} <: AbstractDirichletBCContainer{IT, VT, IV, IM, VV}
-  bookkeeping::BCBookKeeping{IT, IV, IM}
+} <: AbstractBCContainer
+  dofs::IV
+  nodes::IV
   vals::VV
   vals_dot::VV
   vals_dot_dot::VV
@@ -58,7 +58,7 @@ $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
 function DirichletBCContainer(mesh, dof::DofManager, dbc::DirichletBC)
-  bk = BCBookKeeping(mesh, dof, dbc.var_name, dbc.sset_name)
+  bk = BCBookKeeping(mesh, dof, dbc.var_name, sset_name=dbc.sset_name)
 
   # sort nodes and dofs for dirichlet bc
   dof_perm = _unique_sort_perm(bk.dofs)
@@ -72,11 +72,11 @@ function DirichletBCContainer(mesh, dof::DofManager, dbc::DirichletBC)
   vals = zeros(length(bk.nodes))
   vals_dot = zeros(length(bk.nodes))
   vals_dot_dot = zeros(length(bk.nodes))
-  return DirichletBCContainer(bk, vals, vals_dot, vals_dot_dot)
+  return DirichletBCContainer(bk.dofs, bk.nodes, vals, vals_dot, vals_dot_dot)
 end
 
 function Base.length(bc::DirichletBCContainer)
-  return length(bc.bookkeeping.dofs)
+  return length(bc.dofs)
 end
 
 # need checks on if field types are compatable
@@ -87,7 +87,7 @@ based on the stored function
 """
 function _update_bc_values!(bc::DirichletBCContainer, func, X, t, ::KA.CPU)
   ND = num_fields(X)
-  for (n, node) in enumerate(bc.bookkeeping.nodes)
+  for (n, node) in enumerate(bc.nodes)
     X_temp = @views SVector{ND, eltype(X)}(X[:, node])
     bc.vals[n] = func.func(X_temp, t)
     bc.vals_dot[n] = func.func_dot(X_temp, t)
@@ -104,7 +104,7 @@ GPU kernel for updating stored bc values based on the stored function
 KA.@kernel function _update_bc_values_kernel!(bc::DirichletBCContainer, func, X, t)
   I = KA.@index(Global)
   ND = num_fields(X)
-  node = bc.bookkeeping.nodes[I]
+  node = bc.nodes[I]
 
   # hacky for now, but it works
   # can't do X[:, node] on the GPU, this results in a dynamic
@@ -134,22 +134,22 @@ end
 
 # TODO change below names to be more specific to dbcs
 function _update_field_dirichlet_bcs!(U, bc::DirichletBCContainer, ::KA.CPU)
-  for (dof, val) in zip(bc.bookkeeping.dofs, bc.vals)
+  for (dof, val) in zip(bc.dofs, bc.vals)
     U[dof] = val
   end
   return nothing
 end
 
 function _update_field_dirichlet_bcs!(U, V, bc::DirichletBCContainer, ::KA.CPU)
-  for (dof, val, val_dot) in zip(bc.bookkeeping.dofs, bc.vals, bc.vals_dot)
+  for (dof, val, val_dot) in zip(bc.dofs, bc.vals, bc.vals_dot)
     U[dof] = val
     V[dof] = val_dot
   end
   return nothing
 end
 
 function _update_field_dirichlet_bcs!(U, V, A, bc::DirichletBCContainer, ::KA.CPU)
-  for (dof, val, val_dot, val_dot_dot) in zip(bc.bookkeeping.dofs, bc.vals, bc.vals_dot, bc.vals_dot_dot)
+  for (dof, val, val_dot, val_dot_dot) in zip(bc.dofs, bc.vals, bc.vals_dot, bc.vals_dot_dot)
     U[dof] = val
     V[dof] = val_dot
     A[dof] = val_dot_dot
@@ -160,16 +160,15 @@ end
 # COV_EXCL_START
 KA.@kernel function _update_field_dirichlet_bcs_kernel!(U, bc::DirichletBCContainer)
   I = KA.@index(Global)
-  dof = bc.bookkeeping.dofs[I]
-  # val = bc.vals[I]
+  dof = bc.dofs[I]
   U[dof] = bc.vals[I]
 end
 # COV_EXCL_STOP
 
 # COV_EXCL_START
 KA.@kernel function _update_field_dirichlet_bcs_kernel!(U, V, A, bc::DirichletBCContainer)
   I = KA.@index(Global)
-  dof = bc.bookkeeping.dofs[I]
+  dof = bc.dofs[I]
   val = bc.vals[I]
   U[dof] = bc.vals[I]
   V[dof] = bc.vals_dot[I]
@@ -232,7 +231,7 @@ function create_dirichlet_bcs(mesh, dof::DofManager, dirichlet_bcs::Vector{<:Dir
   dirichlet_bcs = DirichletBCContainer.((mesh,), (dof,), dirichlet_bcs)
 
   if length(dirichlet_bcs) > 0
-    temp_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, dirichlet_bcs)
+    temp_dofs = mapreduce(x -> x.dofs, vcat, dirichlet_bcs)
     temp_dofs = unique(sort(temp_dofs))
     update_dofs!(dof, temp_dofs)
   end