Merge pull request #96 from Cthonios/bcs/overhaul

Bcs/overhaul

Author: cmhamel

docs/src/boundary_conditions.md

@@ -1 +1,27 @@
-# Boundary Conditions
+# Boundary Conditions
+This section describes the user facing API for boundary conditions
+along with the implementation details.
+
+## DirichletBC
+We can set up dirichlet boundary conditions on a variable ```u```
+and sideset ```sset_1``` with a zero function as follows.
+```@repl
+using FiniteElementContainers
+bc_func(x, t) = 0.
+bc = DirichletBC(:u, :sset_1, bc_func)
+```
+Internally this is eventually converted in a ```DirichletBCContainer```
+
+### API
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["DirichletBCs.jl"]
+Order = [:type, :function]
+```
+
+## Boundary Condition Implementation Details
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["BoundaryConditions.jl"]
+Order = [:type, :function]
+```

docs/src/index.md

@@ -2,4 +2,25 @@
 CurrentModule = FiniteElementContainers
 ```
 # FiniteElementContainers
-A set of containers to help in the setup of FEM codes.
+```FiniteElementContainers.jl``` is a package whose main purpose is to help 
+researchers develop new finite element method (FEM) applications
+for both well known existing techniques and new novel strategies.
+
+This package is specificially designed with the challenging aspects
+of computational solid mechanics in mind where meshes deform,
+there's path dependence, there's contact between bodies, there are potentially heterogeneous material properties, and other challenges.
+
+If you're primarily interested in writing FEM applications for e.g.
+the Poisson equation or heat equation, there's likely more efficient packages (e.g. ```Gridap``` or ```Ferrite```) out there for this purpose in terms of memory and computational efficiency.
+
+However, if you need to solve problems with multiple material models, meshes where
+there are mixed elements types, etc. this is likely the only julia package
+at the time of writing this that supports such capabilities. 
+
+Inspiration for the software design primarily comes from ```fenics``` and ```MOOSE```.
+We've specifically designed the interface to get around all the shortcomings of ```fenics``` 
+(e.g. boundary conditions are a pain, mixed element types a plain, different blocks are pain
+etc.) 
+
+Our goal is also to ensure all of our methods are next generation hardware capable. This
+means not only supporting things on CPUs but also GPUs (and that doesn't just mean NVIDIA).

ext/FiniteElementContainersAdaptExt.jl

@@ -3,42 +3,37 @@ module FiniteElementContainersAdaptExt
 using Adapt
 using FiniteElementContainers
 
-FiniteElementContainersAdaptExt.cpu(x) = Adapt.adapt_structure(Array, x)
+FiniteElementContainersAdaptExt.cpu(x) = adapt(Array, x)
 
 # Assemblers
 function Adapt.adapt_structure(to, asm::SparseMatrixAssembler)
-  dof = Adapt.adapt_structure(to, asm.dof)
-  pattern = Adapt.adapt_structure(to, asm.pattern)
-  constraint_storage = Adapt.adapt_structure(to, asm.constraint_storage)
-  damping_storage = Adapt.adapt_structure(to, asm.damping_storage)
-  mass_storage = Adapt.adapt_structure(to, asm.mass_storage)
-  residual_storage = Adapt.adapt_structure(to, asm.residual_storage)
-  residual_unknowns = Adapt.adapt_structure(to, asm.residual_unknowns)
-  stiffness_storage = Adapt.adapt_structure(to, asm.stiffness_storage)
   return SparseMatrixAssembler(
-    dof, pattern, 
-    constraint_storage, 
-    damping_storage, mass_storage,
-    residual_storage, residual_unknowns,
-    stiffness_storage
+    adapt(to, asm.dof),
+    adapt(to, asm.pattern),
+    adapt(to, asm.constraint_storage),
+    adapt(to, asm.damping_storage),
+    adapt(to, asm.mass_storage),
+    adapt(to, asm.residual_storage),
+    adapt(to, asm.residual_unknowns),
+    adapt(to, asm.stiffness_storage)
   )
 end
 
 function Adapt.adapt_structure(to, asm::FiniteElementContainers.SparsityPattern)
-  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)
+  Is = adapt(to, asm.Is)
+  Js = adapt(to, asm.Js)
+  unknown_dofs = adapt(to, asm.unknown_dofs)
+  block_sizes = adapt(to, asm.block_sizes)
+  block_offsets = adapt(to, asm.block_offsets)
   #
-  klasttouch = Adapt.adapt_structure(to, asm.klasttouch)
-  csrrowptr = Adapt.adapt_structure(to, asm.csrrowptr)
-  csrcolval = Adapt.adapt_structure(to, asm.csrcolval)
-  csrnzval = Adapt.adapt_structure(to, asm.csrnzval)
+  klasttouch = adapt(to, asm.klasttouch)
+  csrrowptr = adapt(to, asm.csrrowptr)
+  csrcolval = adapt(to, asm.csrcolval)
+  csrnzval = adapt(to, asm.csrnzval)
   #
-  csccolptr = Adapt.adapt_structure(to, asm.csccolptr)
-  cscrowval = Adapt.adapt_structure(to, asm.cscrowval)
-  cscnzval = Adapt.adapt_structure(to, asm.cscnzval)
+  csccolptr = adapt(to, asm.csccolptr)
+  cscrowval = adapt(to, asm.cscrowval)
+  cscnzval = adapt(to, asm.cscnzval)
   return FiniteElementContainers.SparsityPattern(
     Is, Js,
     unknown_dofs, block_sizes, block_offsets,
@@ -48,36 +43,27 @@ function Adapt.adapt_structure(to, asm::FiniteElementContainers.SparsityPattern)
 end
 
 # Boundary Conditions
-function Adapt.adapt_structure(to, bk::FiniteElementContainers.BCBookKeeping{D, S, T, V}) where {D, S, T, V}
-  blocks = Adapt.adapt_structure(to, bk.blocks)
-  dofs = Adapt.adapt_structure(to, bk.dofs)
-  elements = Adapt.adapt_structure(to, bk.elements)
-  nodes = Adapt.adapt_structure(to, bk.nodes)
-  sides = Adapt.adapt_structure(to, bk.sides)
-  return FiniteElementContainers.BCBookKeeping{D, S, T, typeof(blocks)}(blocks, dofs, elements, nodes, sides)
-end
-
-function Adapt.adapt_structure(to, bc::DirichletBC{S, B, F}) where {S, B, F}
-  bk = adapt(to, bc.bookkeeping)
-  func = adapt(to, bc.func)
-  return DirichletBC{S, typeof(bk), typeof(func)}(bk, func)
+function Adapt.adapt_structure(to, bk::FiniteElementContainers.BCBookKeeping{V}) where V
+  blocks = adapt(to, bk.blocks)
+  return FiniteElementContainers.BCBookKeeping{typeof(blocks)}(
+    blocks,
+    adapt(to, bk.dofs),
+    adapt(to, bk.elements),
+    adapt(to, bk.nodes),
+    adapt(to, bk.sides)
+  )
 end
 
 function Adapt.adapt_structure(to, bc::FiniteElementContainers.DirichletBCContainer)
   return FiniteElementContainers.DirichletBCContainer(
     adapt(to, bc.bookkeeping),
-    adapt(to, bc.funcs),
-    adapt(to, bc.func_ids),
+    # adapt(to, bc.funcs),
+    adapt(to, bc.func),
+    # adapt(to, bc.func_ids),
     adapt(to, bc.vals)
   )
 end
 
-# function Adapt.adapt_structure(to, bc::FiniteElementContainers.DirichletBCCollection)
-#   funcs = Adapt.adapt_structure(to, bc.funcs)
-#   func_ids = Adapt.adapt_structure(to, bc.func_ids)
-#   return FiniteElementContainers.DirichletBCCollection(funcs, func_ids)
-# end
-
 # DofManagers
 function Adapt.adapt_structure(to, dof::DofManager{
   T, IDs, 
@@ -88,19 +74,19 @@ function Adapt.adapt_structure(to, dof::DofManager{
   NH1Dofs, NHcurlDofs, NHdivDofs, NL2EDofs, NL2QDofs,
   H1Vars, HcurlVars, HdivVars, L2EVars, L2QVars
 }
-  H1_bc_dofs = Adapt.adapt_structure(to, dof.H1_bc_dofs)
-  H1_unknown_dofs = Adapt.adapt_structure(to, dof.H1_unknown_dofs)
-  Hcurl_bc_dofs = Adapt.adapt_structure(to, dof.Hcurl_bc_dofs)
-  Hcurl_unknown_dofs = Adapt.adapt_structure(to, dof.Hcurl_unknown_dofs)
-  Hdiv_bc_dofs = Adapt.adapt_structure(to, dof.Hdiv_bc_dofs)
-  Hdiv_unknown_dofs = Adapt.adapt_structure(to, dof.Hdiv_unknown_dofs)
-  L2_element_dofs = Adapt.adapt_structure(to, dof.L2_element_dofs)
-  L2_quadrature_dofs = Adapt.adapt_structure(to, dof.L2_quadrature_dofs)
-  H1_vars = Adapt.adapt_structure(to, dof.H1_vars)
-  Hcurl_vars = Adapt.adapt_structure(to, dof.Hcurl_vars)
-  Hdiv_vars = Adapt.adapt_structure(to, dof.Hdiv_vars)
-  L2_element_vars = Adapt.adapt_structure(to, dof.L2_element_vars)
-  L2_quadrature_vars = Adapt.adapt_structure(to, dof.L2_quadrature_vars)
+  H1_bc_dofs = adapt(to, dof.H1_bc_dofs)
+  H1_unknown_dofs = adapt(to, dof.H1_unknown_dofs)
+  Hcurl_bc_dofs = adapt(to, dof.Hcurl_bc_dofs)
+  Hcurl_unknown_dofs = adapt(to, dof.Hcurl_unknown_dofs)
+  Hdiv_bc_dofs = adapt(to, dof.Hdiv_bc_dofs)
+  Hdiv_unknown_dofs = adapt(to, dof.Hdiv_unknown_dofs)
+  L2_element_dofs = adapt(to, dof.L2_element_dofs)
+  L2_quadrature_dofs = adapt(to, dof.L2_quadrature_dofs)
+  H1_vars = adapt(to, dof.H1_vars)
+  Hcurl_vars = adapt(to, dof.Hcurl_vars)
+  Hdiv_vars = adapt(to, dof.Hdiv_vars)
+  L2_element_vars = adapt(to, dof.L2_element_vars)
+  L2_quadrature_vars =adapt(to, dof.L2_quadrature_vars)
   return DofManager{
     T, typeof(H1_bc_dofs), 
     NH1Dofs, NHcurlDofs, NHdivDofs, NL2EDofs, NL2QDofs,
@@ -116,24 +102,24 @@ end
 
 # Fields
 function Adapt.adapt_structure(to, field::L2ElementField{T, NF, V, S}) where {T, NF, V, S}
-  vals = Adapt.adapt_structure(to, field.vals)
+  vals = adapt(to, field.vals)
   return L2ElementField{T, NF, typeof(vals), S}(vals)
 end
 
 function Adapt.adapt_structure(to, field::H1Field{T, NF, V, S}) where {T, NF, V, S}
-  vals = Adapt.adapt_structure(to, field.vals)
+  vals = adapt(to, field.vals)
   return H1Field{T, NF, typeof(vals), S}(vals)
 end
 
 # Function spaces
 function Adapt.adapt_structure(to, fspace::FunctionSpace)
-  coords = Adapt.adapt_structure(to, fspace.coords)
-  elem_conns = Adapt.adapt_structure(to, fspace.elem_conns)
-  elem_id_maps = Adapt.adapt_structure(to, fspace.elem_id_maps)
-  ref_fes = Adapt.adapt_structure(to, fspace.ref_fes)
-  sideset_elems = Adapt.adapt_structure(to, fspace.sideset_elems)
-  sideset_nodes = Adapt.adapt_structure(to, fspace.sideset_nodes)
-  sideset_sides = Adapt.adapt_structure(to, fspace.sideset_sides)
+  coords = adapt(to, fspace.coords)
+  elem_conns = adapt(to, fspace.elem_conns)
+  elem_id_maps = adapt(to, fspace.elem_id_maps)
+  ref_fes = adapt(to, fspace.ref_fes)
+  sideset_elems = adapt(to, fspace.sideset_elems)
+  sideset_nodes = adapt(to, fspace.sideset_nodes)
+  sideset_sides = adapt(to, fspace.sideset_sides)
   return FunctionSpace(
     coords, 
     elem_conns, elem_id_maps, 
@@ -142,78 +128,52 @@ function Adapt.adapt_structure(to, fspace::FunctionSpace)
   )
 end
 
-# # Mesh - won't work, symbols are mutable
-# # function Adapt.adapt_structure(to, mesh::UnstructuredMesh)
-# #   nodal_coords = Adapt.adapt_structure(to, mesh.nodal_coords)
-# #   element_block_names = Adapt.adapt_structure(to, mesh.element_block_names)
-# #   element_types = Adapt.adapt_structure(to, mesh.element_types)
-# #   element_conns = Adapt.adapt_structure(to, mesh.element_conns)
-# #   element_id_maps = Adapt.adapt_structure(to, mesh.element_id_maps)
-# #   nodeset_nodes = Adapt.adapt_structure(to, mesh.nodeset_nodes)
-# #   return UnstructuredMesh(
-# #     nodal_coords, element_block_names,
-# #     element_types, element_conns,
-# #     element_id_maps, nodeset_nodes
-# #   )
-# # end
-
-# # Functions
+# Functions
 function Adapt.adapt_structure(to, var::ScalarFunction)
   syms = names(var)
-  fspace = Adapt.adapt_structure(to, var.fspace)
+  fspace = adapt(to, var.fspace)
   return ScalarFunction{syms, typeof(fspace)}(fspace)
 end
 
 function Adapt.adapt_structure(to, var::SymmetricTensorFunction)
   syms = names(var)
-  fspace = Adapt.adapt_structure(to, var.fspace)
+  fspace = adapt(to, var.fspace)
   return SymmetricTensorFunction{syms, typeof(fspace)}(fspace)
 end
 
 function Adapt.adapt_structure(to, var::TensorFunction)
   syms = names(var)
-  fspace = Adapt.adapt_structure(to, var.fspace)
+  fspace = adapt(to, var.fspace)
   return TensorFunction{syms, typeof(fspace)}(fspace)
 end
 
 function Adapt.adapt_structure(to, var::VectorFunction)
   syms = names(var)
-  fspace = Adapt.adapt_structure(to, var.fspace)
+  fspace = adapt(to, var.fspace)
   return VectorFunction{syms, typeof(fspace)}(fspace)
 end
 
 function Adapt.adapt_structure(to, p::FiniteElementContainers.Parameters)
   return FiniteElementContainers.Parameters(
-    Adapt.adapt_structure(to, p.dirichlet_bcs),
-    Adapt.adapt_structure(to, p.neumann_bcs),
-    Adapt.adapt_structure(to, p.times),
-    Adapt.adapt_structure(to, p.physics),
-    Adapt.adapt_structure(to, p.properties),
-    Adapt.adapt_structure(to, p.state_old),
-    Adapt.adapt_structure(to, p.state_new),
-    Adapt.adapt_structure(to, p.h1_dbcs),
-    Adapt.adapt_structure(to, p.h1_field)
+    adapt(to, p.dirichlet_bcs),
+    adapt(to, p.neumann_bcs),
+    adapt(to, p.times),
+    adapt(to, p.physics),
+    adapt(to, p.properties),
+    adapt(to, p.state_old),
+    adapt(to, p.state_new),
+    adapt(to, p.h1_dbcs),
+    adapt(to, p.h1_field)
   )
 end
 
 function Adapt.adapt_structure(to, p::FiniteElementContainers.TimeStepper)
   return FiniteElementContainers.TimeStepper(
-    Adapt.adapt_structure(to, p.time_start),
-    Adapt.adapt_structure(to, p.time_end),
-    Adapt.adapt_structure(to, p.time_current),
-    Adapt.adapt_structure(to, p.Δt)
+    adapt(to, p.time_start),
+    adapt(to, p.time_end),
+    adapt(to, p.time_current),
+    adapt(to, p.Δt)
   )
 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

ext/FiniteElementContainersCUDAExt.jl

@@ -5,15 +5,17 @@ using CUDA
 using FiniteElementContainers
 using KernelAbstractions
 
-FiniteElementContainers.gpu(x) = Adapt.adapt_structure(CuArray, x)
+FiniteElementContainers.gpu(x) = adapt(CuArray, x)
 
 # this method need to have the assembler initialized first
 # if the stored values in asm.pattern.cscnzval or zero
 # CUDA will error out
 function CUDA.CUSPARSE.CuSparseMatrixCSC(asm::SparseMatrixAssembler)
   @assert typeof(get_backend(asm)) <: CUDABackend "Assembler is not on a CUDA device"
-  @assert length(asm.pattern.cscnzval) > 0 "Need to assemble the assembler once with SparseArrays.sparse!(assembler)"
-  @assert all(x -> x != zero(eltype(asm.pattern.cscnzval)), asm.pattern.cscnzval) "Need to assemble the assembler once with SparseArrays.sparse!(assembler)"
+  # TODO these assert statements are failing for multi dof problems yet
+  # they are doing the right thing
+  # @assert length(asm.pattern.cscnzval) > 0 "Need to assemble the assembler once with SparseArrays.sparse!(assembler)"
+  # @assert all(x -> x != zero(eltype(asm.pattern.cscnzval)), asm.pattern.cscnzval) "Need to assemble the assembler once with SparseArrays.sparse!(assembler)"
   n_dofs = FiniteElementContainers.num_unknowns(asm.dof)
   return CUDA.CUSPARSE.CuSparseMatrixCSC(
     asm.pattern.csccolptr,