neumann bcs now use vectors inste of namedtuples.

Project.toml

@@ -23,8 +23,8 @@ Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
-AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
 AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d"
+AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"

src/Fields.jl

@@ -101,20 +101,30 @@ struct Connectivity{
     nepes::Vector{T}
     nelems::Vector{T}
     offsets::Vector{T}
-end
+# end
+
+    function Connectivity(mats::Vector{<:AbstractArray{<:Integer, 2}})
+        nblocks = length(mats)
+        nepes = map(x -> size(x, 1), mats)
+        nelems = map(x -> size(x, 2), mats)
+        offsets = Vector{eltype(nepes)}(undef, 0)
+        offset = 1
+        for (nepe, nelem) in zip(nepes, nelems)
+            push!(offsets, offset)
+            offset += nepe * nelem
+        end
+        data = mapreduce(vec, vcat, mats)
+        # return Connectivity(data, nblocks, nepes, nelems, offsets)
+        new{eltype(data), typeof(data)}(data, nblocks, nepes, nelems, offsets)
+    end
 
-function Connectivity(mats::Vector{<:AbstractArray{<:Integer, 2}})
-    nblocks = length(mats)
-    nepes = map(x -> size(x, 1), mats)
-    nelems = map(x -> size(x, 2), mats)
-    offsets = Vector{eltype(nepes)}(undef, 0)
-    offset = 1
-    for (nepe, nelem) in zip(nepes, nelems)
-        push!(offsets, offset)
-        offset += nepe * nelem
+    function Connectivity(data, nblocks, nepes, nelems, offsets)
+        new{eltype(data), typeof(data)}(data, nblocks, nepes, nelems, offsets)
+    end
+
+    function Connectivity{T, D}() where {T, D}
+        new{T, D}(T[], 0, T[], T[], T[])
     end
-    data = mapreduce(vec, vcat, mats)
-    return Connectivity(data, nblocks, nepes, nelems, offsets)
 end
 
 function Adapt.adapt_structure(to, conn::Connectivity{T, D}) where {T, D}

src/Parameters.jl

@@ -11,14 +11,17 @@ $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
 struct Parameters{
+  IT       <: Integer,
   RT       <: Number,
-  IV       <: AbstractVector{<:Integer},
+  IV       <: AbstractVector{IT},
   RV       <: AbstractVector{RT},
-  RM       <: AbstractMatrix,
+  RM1      <: AbstractMatrix,
+  RM2      <: AbstractMatrix,
   ICFuncs  <: AbstractVector,
   DBCFuncs <: AbstractVector,
   SrcFuncs <: AbstractVector,
-  NBCs,
+  # NBCs,
+  NBCFuncs <: AbstractVector,
   RBCs,
   Phys,
   Props,
@@ -27,9 +30,10 @@ struct Parameters{
 } <: AbstractParameters
   ics::InitialConditions{ICFuncs, IV, RV}
   dirichlet_bcs::DirichletBCs{DBCFuncs, IV, RV}
-  neumann_bcs::NBCs
+  # neumann_bcs::NBCs
+  neumann_bcs::NeumannBCs{NBCFuncs, IT, IV, RM1}
   robin_bcs::RBCs
-  sources::Sources{SrcFuncs, RM}
+  sources::Sources{SrcFuncs, RM2}
   times::TimeStepper{RT}
   physics::Phys
   properties::Props
@@ -241,7 +245,7 @@ function update_bc_values!(p::AbstractParameters, assembler)
   X = coordinates(p)
   t = current_time(p)
   update_bc_values!(p.dirichlet_bcs, X, t)
-  update_bc_values!(p.neumann_bcs, X, t)
+  update_bc_values!(p.neumann_bcs, assembler, X, t)
   update_source_values!(p.sources, assembler, X, t)
 
   # TODO how to handle Robin BCs?

src/assemblers/WeaklyEnforcedBCs.jl

@@ -8,7 +8,7 @@ function assemble_vector_neumann_bc!(
   # been conducted previously?
   _update_for_assembly!(p, assembler.dof, Uu)
   assemble_vector_weakly_enforced_bc!(
-    assembler.residual_storage,
+    assembler.residual_storage, assembler.dof,
     p.field, coordinates(p), p.neumann_bcs
   )
   return nothing
@@ -25,7 +25,7 @@ function assemble_vector_robin_bc!(
   _update_for_assembly!(p, assembler.dof, Uu)
   update_bc_values!(p.robin_bcs, coordinates(p), current_time(p), p.field)
   assemble_vector_weakly_enforced_bc!(
-    assembler.residual_storage,
+    assembler.residual_storage, assembler.dof,
     p.field, coordinates(p), p.neumann_bcs
   )
   return nothing
@@ -36,14 +36,16 @@ end
 $(TYPEDSIGNATURES)
 """
 function assemble_vector_weakly_enforced_bc!(
-  storage, U, X, bcs
+  storage, dof, U, X, bcs
 )
   # do not zero!
-  for bc in values(bcs.bc_caches)
+  for (n, bc) in enumerate(values(bcs.bc_caches))
+    block_id = bcs.block_ids[n]
+    ref_fe = block_reference_element(function_space(dof), block_id)
     _assemble_block_vector_weakly_enforced_bc!(
       storage,
       U, X,
-      bc.element_conns.data, bc.ref_fe, bc.sides, bc.vals
+      bc.element_conns.data, ref_fe, bc.sides, bc.vals
     )
   end
 end

src/bcs/BoundaryConditions.jl

@@ -175,20 +175,9 @@ $(TYPEDFIELDS)
 abstract type AbstractWeaklyEnforcedBCContainer{
   IT <: Integer,
   IV <: AbstractArray{IT, 1},
-  RV <: AbstractArray{<:Union{<:Number, <:SVector}, 2},
-  RE <: ReferenceFE
+  RV <: AbstractArray{<:Union{<:Number, <:SVector}, 2}
 } <: AbstractBCContainer{IV, RV} end
 
-function Adapt.adapt_structure(to, bc::AbstractWeaklyEnforcedBCContainer)
-  el_conns = adapt(to, bc.element_conns)
-  elements = adapt(to, bc.elements)
-  sides = adapt(to, bc.sides)
-  ref_fe = adapt(to, bc.ref_fe)
-  vals = adapt(to, bc.vals)
-  type = eval(typeof(bc).name.name)
-  return type(el_conns, elements, sides, ref_fe, vals)
-end
-
 Base.length(bc::AbstractWeaklyEnforcedBCContainer) = size(bc.vals, 2)
 
 function Base.show(io::IO, bc::AbstractWeaklyEnforcedBCContainer)
@@ -231,15 +220,6 @@ end
 """
 function update_bc_values! end
 
-
-function Adapt.adapt_structure(to, bcs::AbstractBCs)
-  type = typeof(bcs).name.name
-  return eval(type)(
-    map(x -> adapt(to, x), bcs.bc_caches),
-    adapt(to, bcs.bc_funcs)
-  )
-end
-
 Base.length(bcs::AbstractBCs) = length(bcs.bc_funcs)
 
 # function Base.show(io::IO, bcs::AbstractBCs)
@@ -273,6 +253,10 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type)
   fspace = function_space(dof)
   new_bcs = type[]
   new_funcs = Function[]
+  block_ids = Int[]
+  block_names = String[]
+  # sideset_ids = Int[]
+  # sideset_names = String[]
   for (bk, func) in zip(bks, funcs)
     blocks = sort(unique(bk.blocks))
 
@@ -283,6 +267,10 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type)
 
     for block in blocks
       block_name = mesh.element_block_names_map[block]
+      block_id = findfirst(x -> x == block_name, mesh.element_block_names)
+      push!(block_ids, block_id)
+      push!(block_names, block_name)
+
       ids = findall(x -> x == block, bk.blocks)
       new_blocks = bk.blocks[ids]
       new_elements = bk.elements[ids]
@@ -311,13 +299,13 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type)
 
       if type <: NeumannBCContainer
         new_bc = type(
-          conns, new_bk.elements, new_bk.sides, ref_fe, vals
+          conns, new_bk.elements, new_bk.sides, vals
         )
       elseif type <: RobinBCContainer
         # dvalsdu = copy(vals)
         dvalsdu = zeros(SMatrix{ND, ND, Float64, ND * ND}, NQ, length(bk.sides))
         new_bc = type(
-          conns, new_bk.elements, new_bk.sides, ref_fe, vals, dvalsdu
+          conns, new_bk.elements, new_bk.sides, vals, dvalsdu
         )
       else
         _unsupported_weak_bc_error("Unsupported weak bc type $type encountered in _setup_weakly_enforced_bc_container")
@@ -326,5 +314,5 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type)
       push!(new_funcs, func)
     end
   end
-  return new_bcs, new_funcs
+  return new_bcs, new_funcs, block_ids, block_names
 end

src/bcs/DirichletBCs.jl

@@ -338,7 +338,20 @@ struct DirichletBCs{
 end
 
 function Adapt.adapt_structure(to, bcs::DirichletBCs{BCF, IV, RV}) where {BCF, IV, RV}
-  bc_cache = adapt(to, bcs.bc_cache)
+  # NOTE
+  # below logic is needed due to improper
+  # adapt mapping for an empty array in julia 1.10/1.11
+  # where Vector{T}(undef, 0) gets mappend to Vector{Any}
+  if length(bcs.bc_lengths) > 0
+    bc_cache = adapt(to, bcs.bc_cache)
+  else
+    temp_int = adapt(to, zeros(Int, 0))
+    temp_floats = adapt(to, zeros(Float64, 0))
+    bc_cache = DirichletBCContainer{typeof(temp_int), typeof(temp_float)}(
+      temp_int, copy(temp_int),
+      temp_floats, copy(temp_floats), copy(temp_floats)
+    )
+  end
   return DirichletBCs{BCF, typeof(bc_cache.dofs), typeof(bc_cache.vals)}(
     bc_cache,
     bcs.bc_funcs,

src/bcs/NeumannBCs.jl

@@ -28,14 +28,28 @@ Internal implementation of dirichlet BCs
 struct NeumannBCContainer{
   IT <: Integer,
   IV <: AbstractArray{IT, 1},
-  RV <: AbstractArray{<:Union{<:Number, <:SVector}, 2},
-  RE <: ReferenceFE
-} <: AbstractWeaklyEnforcedBCContainer{IT, IV, RV, RE}
+  RM <: AbstractArray{<:Union{<:Number, <:SVector}, 2}
+} <: AbstractWeaklyEnforcedBCContainer{IT, IV, RM}
   element_conns::Connectivity{IT, IV}
   elements::IV
   sides::IV
-  ref_fe::RE
-  vals::RV
+  vals::RM
+
+  function NeumannBCContainer(element_conns::Connectivity{IT, IV}, elements::IV, sides::IV, vals::RM) where {IT, IV, RM}
+    new{IT, IV, RM}(element_conns, elements, sides, vals)
+  end
+
+  function NeumannBCContainer{IT, IV, RM}() where {IT, IV, RM}
+    new{IT, IV, RM}(Connectivity{IT, IV}(), IV(undef, 0), IV(undef, 0), RM(undef, 0))
+  end
+end
+
+function Adapt.adapt_structure(to, bc::NeumannBCContainer)
+  el_conns = adapt(to, bc.element_conns)
+  elements = adapt(to, bc.elements)
+  sides = adapt(to, bc.sides)
+  vals = adapt(to, bc.vals)
+  return NeumannBCContainer(el_conns, elements, sides, vals)
 end
 
 function _update_bc_values!(vals, func, ref_fe, conns, sides, X, t)
@@ -51,12 +65,24 @@ function _update_bc_values!(vals, func, ref_fe, conns, sides, X, t)
   end
 end
 
+struct NeumannBCFunction{F} <: AbstractBCFunction{F}
+  func::F
+end
+
+function (f::NeumannBCFunction)(x, t)
+  return f.func(x, t)
+end
+
 struct NeumannBCs{
-  BCCaches <: NamedTuple, 
-  BCFuncs  <: NamedTuple
+  BCFuncs,
+  IT      <: Integer,
+  IV      <: AbstractArray{IT, 1},
+  RM      <: AbstractArray{<:Union{<:Number, <:SVector}, 2},
 } <: AbstractBCs{BCFuncs}
-  bc_caches::BCCaches
+  bc_caches::Vector{NeumannBCContainer{IT, IV, RM}}
   bc_funcs::BCFuncs
+  block_ids::Vector{Int}
+  block_names::Vector{String}
 end
 
 # note this method has the potential to make 
@@ -69,21 +95,50 @@ end
 # TODO below method also currently likely doesn't
 # handle blocks correclty 
 function NeumannBCs(mesh, dof::DofManager, neumann_bcs::Vector{NeumannBC})
+  bc_caches = NeumannBCContainer{Int, Vector{Int}, Matrix{Float64}}[]
   if length(neumann_bcs) == 0
-    return NeumannBCs(NamedTuple(), NamedTuple())
+    return NeumannBCs(bc_caches, NeumannBCFunction[], Int[], String[])
   end
 
-  new_bcs, new_funcs = _setup_weakly_enforced_bc_container(mesh, dof, neumann_bcs, NeumannBCContainer)
+  new_bcs, new_funcs, block_ids, block_names = _setup_weakly_enforced_bc_container(mesh, dof, neumann_bcs, NeumannBCContainer)
 
-  syms = tuple(map(x -> Symbol("neumann_bc_$x"), 1:length(new_bcs))...)
-  new_bcs = NamedTuple{syms}(tuple(new_bcs...))
-  new_funcs = NamedTuple{syms}(tuple(new_funcs...))
-  return NeumannBCs(new_bcs, new_funcs)
+  # TODO fix this up eventually
+  new_bcs = convert(Vector{typeof(new_bcs[1])}, new_bcs)
+  new_funcs = NeumannBCFunction.(new_funcs)
+  # syms = tuple(map(x -> Symbol("neumann_bc_$x"), 1:length(new_bcs))...)
+  # new_bcs = NamedTuple{syms}(tuple(new_bcs...))
+  # new_funcs = NamedTuple{syms}(tuple(new_funcs...))
+  return NeumannBCs(new_bcs, new_funcs, block_ids, block_names)
+end
+
+function Adapt.adapt_structure(to, bcs::NeumannBCs)
+  # NOTE
+  # below logic is needed due to improper
+  # adapt mapping for an empty array in julia 1.10/1.11
+  # where Vector{T}(undef, 0) gets mappend to Vector{Any}
+  if length(bcs.bc_caches) > 0
+    bc_caches = map(x -> adapt(to, x), bcs.bc_caches)
+  else
+    temp_int = adapt(to, zeros(Int, 0))
+    temp_floats = adapt(to, zeros(Float64, 0, 0))
+    bc_caches = NeumannBCContainer{Int, typeof(temp_int), typeof(temp_floats)}[]
+  end
+  return NeumannBCs(
+    bc_caches,
+    bcs.bc_funcs,
+    bcs.block_ids,
+    bcs.block_names
+  )
 end
 
-function update_bc_values!(bcs::NeumannBCs, X, t)
-  for (bc, func) in zip(bcs.bc_caches, bcs.bc_funcs)
-    _update_bc_values!(bc.vals, func, bc.ref_fe, bc.element_conns.data, bc.sides, X, t)
+function update_bc_values!(bcs::NeumannBCs, assembler, X, t)
+  fspace = function_space(assembler.dof)
+  for n in axes(bcs.block_ids, 1)
+    bc = bcs.bc_caches[n]
+    block_id = bcs.block_ids[n]
+    func = bcs.bc_funcs[n]
+    ref_fe = block_reference_element(fspace, block_id)
+    _update_bc_values!(bc.vals, func, ref_fe, bc.element_conns.data, bc.sides, X, t)
   end
   return nothing
 end