Merge pull request #172 from Cthonios/bcs/dirichlet-rework

reworking dirichlet bcs to hold onto velocity and acceleration values…

Author: cmhamel

Project.toml

@@ -1,11 +1,12 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
 authors = ["Craig M. Hamel <[email protected]> and contributors"]
-version = "0.8.1"
+version = "0.9.0"
 
 [deps]
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -57,9 +58,8 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 
 [targets]
-test = ["Adapt", "AMDGPU", "Aqua", "CUDA", "Exodus", "ForwardDiff", "Test", "TestSetExtensions"]
+test = ["Adapt", "AMDGPU", "Aqua", "CUDA", "Exodus", "Test", "TestSetExtensions"]

ext/FiniteElementContainersAdaptExt.jl

@@ -49,26 +49,22 @@ end
 
 # Boundary Conditions
 function Adapt.adapt_structure(to, bk::FiniteElementContainers.BCBookKeeping{V}) where V
-  blocks = adapt(to, bk.blocks)
-  side_nodes = adapt(to, bk.side_nodes)
-  return FiniteElementContainers.BCBookKeeping{typeof(side_nodes), typeof(blocks)}(
-    blocks,
+  return FiniteElementContainers.BCBookKeeping(
+    adapt(to, bk.blocks),
     adapt(to, bk.dofs),
     adapt(to, bk.elements),
     adapt(to, bk.nodes),
     adapt(to, bk.sides),
-    side_nodes
+    adapt(to, bk.side_nodes)
   )
 end
 
 function Adapt.adapt_structure(to, bc::FiniteElementContainers.DirichletBCContainer)
   bk = adapt(to, bc.bookkeeping)
   vals = adapt(to, bc.vals)
-  return FiniteElementContainers.DirichletBCContainer{
-    typeof(bk), eltype(vals), typeof(vals)
-  }(
-    bk, 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)
 end
 
 function Adapt.adapt_structure(to, bc::FiniteElementContainers.NeumannBCContainer)
@@ -77,12 +73,7 @@ function Adapt.adapt_structure(to, bc::FiniteElementContainers.NeumannBCContaine
   surf_conns = adapt(to, bc.surface_conns)
   ref_fe = adapt(to, bc.ref_fe)
   vals = adapt(to, bc.vals)
-  return FiniteElementContainers.NeumannBCContainer{
-    typeof(bk), typeof(el_conns), typeof(surf_conns),
-    typeof(ref_fe), eltype(vals), typeof(vals)
-  }(
-    bk, el_conns, surf_conns, ref_fe, vals
-  )
+  return FiniteElementContainers.NeumannBCContainer(bk, el_conns, surf_conns, ref_fe, vals)
 end
 
 # DofManagers
@@ -116,11 +107,7 @@ function Adapt.adapt_structure(to, fspace::FunctionSpace)
   coords = adapt(to, fspace.coords)
   elem_conns = adapt(to, fspace.elem_conns)
   ref_fes = adapt(to, fspace.ref_fes)
-  return FunctionSpace(
-    coords, 
-    elem_conns, 
-    ref_fes
-  )
+  return FunctionSpace(coords, elem_conns, ref_fes)
 end
 
 function Adapt.adapt_structure(to, var::T) where T <: FiniteElementContainers.AbstractFunction

src/FiniteElementContainers.jl

@@ -32,8 +32,8 @@ export L2QuadratureField
 export DofManager
 export create_field
 export create_unknowns
-export update_bcs!
 export update_dofs!
+export update_field_dirichlet_bcs!
 export update_field_unknowns!
 
 # Formulations
@@ -133,6 +133,7 @@ export MappedSurfaceInterpolants
 import KernelAbstractions as KA
 using Atomix
 using DocStringExtensions
+using ForwardDiff
 using Krylov
 using LinearAlgebra
 using ReferenceFiniteElements

src/Parameters.jl

@@ -113,7 +113,6 @@ function Parameters(
   dirichlet_bcs, dirichlet_bc_funcs = create_dirichlet_bcs(
     mesh, assembler.dof, dirichlet_bcs
   )
-  # update_dofs!(assembler.dof, )
 
   if neumann_bcs !== NamedTuple()
     syms = map(x -> Symbol("neumann_bc_$x"), 1:length(neumann_bcs))
@@ -184,13 +183,6 @@ function create_parameters(
   return Parameters(mesh, assembler, physics, props, dirichlet_bcs, neumann_bcs, times)
 end
 
-function update_bcs!(p::Parameters)
-  for bc in values(p.dirichlet_bcs)
-    _update_bcs!(bc, p.h1_field, KA.get_backend(bc))
-  end
-  return nothing
-end
-
 """
 $(TYPEDSIGNATURES)
 This method is used to update the stored bc values.
@@ -213,13 +205,13 @@ function update_dofs!(asm::SparseMatrixAssembler, p::Parameters)
 end
 
 function _update_for_assembly!(p::Parameters, dof::DofManager, Uu)
-  update_bcs!(p)
+  update_field_dirichlet_bcs!(p.h1_field, p.dirichlet_bcs)
   update_field_unknowns!(p.h1_field, dof, Uu)
   return nothing
 end
 
 function _update_for_assembly!(p::Parameters, dof::DofManager, Uu, Vu)
-  update_bcs!(p)
+  update_field_dirichlet_bcs!(p.h1_field, p.dirichlet_bcs)
   update_field_unknowns!(p.h1_field, dof, Uu)
   update_field_unknowns!(p.h1_hvp_scratch_field, dof, Vu)
   return nothing

src/assemblers/QuadratureQuantity.jl

@@ -20,8 +20,7 @@ function assemble_quadrature_quantity!(
   fspace = function_space(dof)
   t = current_time(p.times)
   Δt = time_step(p.times)
-  update_bcs!(p)
-  update_field_unknowns!(p.h1_field, dof, Uu)
+  _update_for_assembly!(p, dof, Uu)
   for (b, (field, conns, block_physics, state_old, state_new, props)) in enumerate(zip(
     values(storage),
     values(fspace.elem_conns), 

src/assemblers/SparseMatrixAssembler.jl

@@ -115,26 +115,16 @@ $(TYPEDSIGNATURES)
 TODO add symbol to interface
 """
 function SparseArrays.sparse!(
-  assembler::SparseMatrixAssembler, ::Val{:stiffness};
-  unknowns_only=true
+  assembler::SparseMatrixAssembler, ::Val{:stiffness}
 )
   pattern = assembler.pattern
   storage = assembler.stiffness_storage
-  if unknowns_only
-    return @views SparseArrays.sparse!(
-      pattern.Is, pattern.Js, storage[assembler.pattern.unknown_dofs],
-      length(pattern.klasttouch), length(pattern.klasttouch), +, pattern.klasttouch,
-      pattern.csrrowptr, pattern.csrcolval, pattern.csrnzval,
-      pattern.csccolptr, pattern.cscrowval, pattern.cscnzval
-    )
-  # else
-  #   return SparseArrays.sparse!(
-  #     pattern.Is, pattern.Js, storage[assembler.pattern.unknown_dofs],
-  #     length(pattern.klasttouch), length(pattern.klasttouch), +, pattern.klasttouch,
-  #     pattern.csrrowptr, pattern.csrcolval, pattern.csrnzval,
-  #     pattern.csccolptr, pattern.cscrowval, pattern.cscnzval
-  #   )
-  end
+  return @views SparseArrays.sparse!(
+    pattern.Is, pattern.Js, storage[assembler.pattern.unknown_dofs],
+    length(pattern.klasttouch), length(pattern.klasttouch), +, pattern.klasttouch,
+    pattern.csrrowptr, pattern.csrcolval, pattern.csrnzval,
+    pattern.csccolptr, pattern.cscrowval, pattern.cscnzval
+  )
 end
 
 """
@@ -152,13 +142,8 @@ function SparseArrays.sparse!(assembler::SparseMatrixAssembler, ::Val{:mass})
   )
 end
 
-function SparseArrays.spdiagm(assembler::SparseMatrixAssembler)
-  return SparseArrays.spdiagm(assembler.constraint_storage)
-end
-
 function constraint_matrix(assembler::SparseMatrixAssembler)
-  # TODO specialize to CPU/GPU
-  return SparseArrays.spdiagm(assembler)
+  return Diagonal(assembler.constraint_storage)
 end
 
 function _mass(assembler::SparseMatrixAssembler, ::KA.CPU)

src/bcs/BoundaryConditions.jl

@@ -1,3 +1,26 @@
+struct VariableNameNotFoundError <: Exception
+end
+
+_var_not_found_err() = throw(VariableNameNotFoundError())
+
+function _dof_index_from_var_name(dof, var_name)
+  dof_index = 0
+  found = false
+  for name in names(dof.var)
+      dof_index = dof_index + 1
+      if var_name == name
+          found = true
+          break
+      end
+  end
+
+  if dof_index == 0 || dof_index > length(names(dof.var))
+    _var_not_found_err()
+  end
+
+  return dof_index
+end
+
 function _unique_sort_perm(array::AbstractArray{T, 1}) where T <: Number
   # chatgpt BS below. Make this not use a dict
   sorted_unique = sort(unique(array))
@@ -15,8 +38,9 @@ for all types of boundary conditions.
 TODO need to add a domain ID for extending to Schwarz
 """
 struct BCBookKeeping{
-  M <: AbstractArray{<:Integer, 2},
-  V <: AbstractArray{<:Integer, 1}
+  I <: Integer,
+  V <: AbstractArray{I, 1},
+  M <: AbstractArray{I, 2}
 }
   blocks::V
   dofs::V
@@ -33,36 +57,9 @@ $(TYPEDSIGNATURES)
 function BCBookKeeping(
   mesh, dof::DofManager, var_name::Symbol, sset_name::Symbol
 )
-  # check if var exists
-  # var_index = 0
-  dof_index = 0
-  found = false
-  # for (vi, var) in enumerate(dof.H1_vars)
-  #   for name in names(var)
-  #     dof_index = dof_index + 1
-  #     if var_name == name
-  #       var_index = vi
-  #       found = true
-  #       break
-  #     end
-  #   end
-  # end
-  for name in names(dof.var)
-    dof_index = dof_index + 1
-    if var_name == name
-      # var_index = vi
-      found = true
-      break
-    end
-  end
-
-  # some error checking
-  # @assert found == true "Failed to find variable $var_name"
-  # @assert dof_index <= length(mapreduce(x -> names(x), +, dof.H1_vars)) "Found invalid dof index"
-  @assert dof_index <= length(names(dof.var)) "Found invalid dof index"
+  # get dof index associated with this var
+  dof_index = _dof_index_from_var_name(dof, var_name)
 
-  # TODO
-  # fspace = dof.H1_vars[var_index].fspace
   fspace = function_space(dof)
 
   # get sset specific fields
@@ -71,7 +68,6 @@ function BCBookKeeping(
   sides = getproperty(mesh.sideset_sides, sset_name)
   side_nodes = getproperty(mesh.sideset_side_nodes, sset_name)
 
-  dofs = Vector{Int64}(undef, 0)
   blocks = Vector{Int64}(undef, 0)
 
   # gather the blocks that are present in this sideset
@@ -86,48 +82,11 @@ function BCBookKeeping(
     end
   end
 
-  # setting up dofs for use in dirichlet bcs
-  # TODO only works for H1
-  # ND, NN = num_dofs_per_node(dof), num_nodes(dof)
-  ND, NN = size(dof)
-  n_total_dofs = ND * NN
-  all_dofs = reshape(1:n_total_dofs, ND, NN)
-  # temp_dofs = all_dofs[dof_index, temp_nodes]
-  temp_dofs = all_dofs[dof_index, nodes]
-  append!(dofs, temp_dofs)
-  # end
-
-  # now sort and unique this stuff
-  # dof_perm = _unique_sort_perm(dofs)
-  # # el_perm = _unique_sort_perm(elements)
-
-  # # do permutations
-  # # TODO need to do this in different BC types
-  # # since different BC types might want things
-  # # organized differently
-  # # blocks_new = blocks[el_perm]
-  # dofs_new = dofs[dof_perm]
-  # # elements_new = elements[el_perm]
-  # nodes_new = nodes[dof_perm]
-  # # sides_new = sides[el_perm]
-  # # side_nodes_new = side_nodes[:, el_perm]
-
-  # # resize!(blocks, length(blocks_new))
-  # resize!(dofs, length(dofs_new))
-  # # resize!(elements, length(elements_new))
-  # resize!(nodes, length(nodes_new))
-  # # resize!(sides, length(sides_new))
-  # # resize!(side_nodes, size(side_nodes_new)...)
-  # # reshape(side_nodes, size(side_nodes_new))
-
-  # # copyto!(blocks, blocks_new)
-  # copyto!(dofs, dofs_new)
-  # # copyto!(elements, elements_new)
-  # copyto!(nodes, nodes_new)
-  # # copyto!(sides, sides_new)
-  # # copyto!(side_nodes, side_nodes_new)
-
-  return BCBookKeeping{typeof(side_nodes), typeof(blocks)}(
+  # setup dofs local to this BC
+  all_dofs = reshape(1:length(dof), size(dof))
+  dofs = all_dofs[dof_index, nodes]
+
+  return BCBookKeeping(
     blocks, dofs, elements, nodes, sides, side_nodes
   )
 end
@@ -158,10 +117,12 @@ $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
 abstract type AbstractBCContainer{
-  B <: BCBookKeeping,
-  T <: Union{<:Number, <:SVector},
+  IT <: Integer,
+  VT <: Union{<:Number, <:SVector},
   N,
-  V <: AbstractArray{T, N}
+  IV <: AbstractArray{IT, 1},
+  IM <: AbstractArray{IT, 2},
+  VV <: AbstractArray{VT, N}
 } end
 
 KA.get_backend(x::AbstractBCContainer) = KA.get_backend(x.vals)
@@ -171,6 +132,13 @@ function Base.show(io::IO, bc::AbstractBCContainer)
   println(io, "$(bc.bookkeeping)")
 end
 
+"""
+$(TYPEDEF)
+$(TYPEDSIGNATURES)
+$(TYPEDFIELDS)
+"""
+abstract type AbstractBCFunction{F} end
+
 """
 $(TYPEDSIGNATURES)
 Wrapper that is generic for all architectures to
@@ -181,8 +149,6 @@ function update_bc_values!(bcs, funcs, X, t)
     backend = KA.get_backend(bc)
     _update_bc_values!(bc, func, X, t, backend)
   end
-  backend = KA.get_backend(X)
-  KA.synchronize(backend)
   return nothing
 end