Merge pull request #164 from Cthonios/fspace-cleanup

removing sideset bookkeeping stuff from function spaces since this is…

Author: cmhamel

ext/FiniteElementContainersAdaptExt.jl

@@ -113,19 +113,11 @@ end
 function Adapt.adapt_structure(to, fspace::FunctionSpace)
   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)
-  sideset_side_nodes = adapt(to, fspace.sideset_side_nodes)
   return FunctionSpace(
     coords, 
-    elem_conns, elem_id_maps, 
-    ref_fes,
-    sideset_elems, 
-    sideset_nodes, 
-    sideset_sides, sideset_side_nodes
+    elem_conns, 
+    ref_fes
   )
 end
 

src/FiniteElementContainers.jl

@@ -83,7 +83,6 @@ export num_dimensions
 export num_fields
 export num_nodes
 export num_nodes_per_element
-export num_q_points
 export sideset
 export sideset_ids
 export sideset_names

src/FunctionSpaces.jl

@@ -13,18 +13,12 @@ $(TYPEDFIELDS)
 """
 struct FunctionSpace{
   Coords,
-  ElemConns, ElemIdMaps, 
-  RefFEs,
-  SSetIVs, SSetIMs
+  ElemConns,
+  RefFEs
 } <: AbstractFunctionSpace
   coords::Coords
   elem_conns::ElemConns
-  elem_id_maps::ElemIdMaps
   ref_fes::RefFEs
-  sideset_elems::SSetIVs
-  sideset_nodes::SSetIVs
-  sideset_sides::SSetIVs
-  sideset_side_nodes::SSetIMs
 end
 
 function _setup_ref_fes(mesh::AbstractMesh, interp_type, q_degree)
@@ -69,11 +63,8 @@ function FunctionSpace(mesh::AbstractMesh, ::Type{H1Field}, interp_type, q_degre
 
   return FunctionSpace(
     coords, 
-    mesh.element_conns, mesh.element_id_maps, 
-    ref_fes,
-    mesh.sideset_elems, 
-    mesh.sideset_nodes, 
-    mesh.sideset_sides, mesh.sideset_side_nodes
+    mesh.element_conns,
+    ref_fes
   )
 end
 
@@ -100,11 +91,8 @@ function FunctionSpace(mesh::AbstractMesh, ::Type{L2ElementField}, interp_type,
 
   return FunctionSpace(
     coords,
-    nothing, nothing, # TODO what makes sense here?
-    ref_fes,
-    mesh.sideset_elems, 
-    mesh.sideset_nodes, 
-    mesh.sideset_sides, mesh.sideset_side_nodes
+    nothing, # TODO what makes sense here?
+    ref_fes
   )
 end
 
@@ -123,24 +111,13 @@ function FunctionSpace(mesh::AbstractMesh, ::Type{L2QuadratureField}, interp_typ
     push!(coords_vals, coords_temp)
   end
 
-  if num_fields(mesh.nodal_coords) == 1
-    temp_syms = (:coords_x,)
-  elseif num_fields(mesh.nodal_coords) == 2
-    temp_syms = (:coords_x, :coords_y)
-  else
-    temp_syms = (:coords_x, :coords_y, :coords_z)
-  end
-
   coords_vals = L2QuadratureField.(coords_vals)
   coords = NamedTuple{tuple(coords_syms...)}(tuple(coords_vals...))
 
   return FunctionSpace(
     coords,
-    nothing, nothing, # TODO what makes sense here?
-    ref_fes,
-    mesh.sideset_elems, 
-    mesh.sideset_nodes, 
-    mesh.sideset_sides, mesh.sideset_side_nodes
+    nothing, # TODO what makes sense here?
+    ref_fes
   )
 end
 
@@ -153,68 +130,12 @@ function Base.show(io::IO, fspace::FunctionSpace)
   println(io, "  Type: $(typeof(fspace.coords).name.name)")
   for (key, ref_fe) in enumerate(fspace.ref_fes)
     println(io, "    Block: $key")
-    # println(io, "      Number of elements: $(fspace)")
-    # println(io, "  $ref_fe")
   end
 end
 
-function connectivity(fspace, block)
-  return Base.getproperty(fspace.elem_conns, block)
-end
-
-function connectivity(fspace, e::Int, block::Symbol)
-  temp = @views Base.getproperty(fspace.elem_conns, block)[:, e]
-  # NN = num_vertices(fspace.ref_fes[block])
-  # return SVector{NN, eltype(temp)}(temp)
-  return temp
-end
-
-function connectivity(fspace, e::Int, block_num::Int)
-  block_sym = keys(fspace.elem_conns)[block_num]
-  return @views Base.getproperty(fspace.elem_conns, block_sym)[:, e]
-end
-
-# connectivity(fspace, e, block::Symbol) = _connectivity(fspace, e, Val(block))
-
-function connectivity(fspace, n::Int, e::Int, block)
-  return @views Base.getproperty(fspace.elem_conns, block)[n, e]
-end 
-
-coordinates(fspace::FunctionSpace) = fspace.coords
-
-# function coordinates(fspace::FunctionSpace, e, block)
-#   conn = _connectivity(fspace, e, block)
-#   return @views fspace.coords[:, conn]
-# end
-
-function dof_connectivity(fspace, e, block, n_dofs)
-  ids = reshape(1:length(fspace.coords), size(fspace.coords)...)
-  block_name = keys(fspace.elem_conns)[block]
-  conns = getproperty(fspace.elem_conns, block_name)
-  # dof_conn = @views reshape(ids[:, conns[:, e]], n_dofs * size(conns, 1))
-  dof_conn = reshape(ids[:, conns[:, e]], n_dofs * size(conns, 1))
-  return dof_conn
-end
-
 function map_shape_function_gradients(X, ∇N_ξ)
   J     = (X * ∇N_ξ)'
   J_inv = inv(J)
   ∇N_X  = (J_inv * ∇N_ξ')'
   return ∇N_X
 end
-
-function num_dimensions(fspace::FunctionSpace, block_sym)
-  return ReferenceFiniteElements.dimension(getproperty(fspace.ref_fes, block_sym))
-end
-
-function num_elements(fspace::FunctionSpace, block_sym)
-  return num_elements(getproperty(fspace.elem_conns, block_sym))
-end
-
-function num_nodes_per_element(fspace::FunctionSpace, block_sym)
-  return ReferenceFiniteElements.num_vertices(getproperty(fspace.ref_fes, block_sym))
-end 
-
-function num_q_points(fspace::FunctionSpace, block_sym)
-  return ReferenceFiniteElements.num_quadrature_points(getproperty(fspace.ref_fes, block_sym))
-end

src/Parameters.jl

@@ -40,7 +40,7 @@ end
 # 4. add different fspace types
 # 5. convert vectors of dbcs/nbcs into namedtuples - done
 function Parameters(
-  assembler, physics,
+  mesh, assembler, physics,
   properties,
   dirichlet_bcs, 
   neumann_bcs, 
@@ -86,11 +86,9 @@ function Parameters(
     # create state variables for this block physics
     NS = num_states(val)
     NQ = ReferenceFiniteElements.num_quadrature_points(
-      # getfield(values(assembler.dof.H1_vars)[1].fspace.ref_fes, key)
       getfield(function_space(assembler.dof).ref_fes, key)
     )
     NE = size(
-      # getfield(values(assembler.dof.H1_vars)[1].fspace.elem_conns, key),
       getfield(function_space(assembler.dof).elem_conns, key),
       2
     )
@@ -108,20 +106,20 @@ function Parameters(
 
   # properties = NamedTuple{keys(physics)}(tuple(properties...))
 
-  state_new = copy(state_old)
+  state_new = deepcopy(state_old)
   state_old = NamedTuple{keys(physics)}(tuple(state_old...))
   state_new = NamedTuple{keys(physics)}(tuple(state_new...))
 
   dirichlet_bcs, dirichlet_bc_funcs = create_dirichlet_bcs(
-    assembler.dof, 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))
 
     if length(neumann_bcs) > 1
-      neumann_bcs, neumann_bc_funcs = create_neumann_bcs(assembler.dof, neumann_bcs)
+      neumann_bcs, neumann_bc_funcs = create_neumann_bcs(mesh, assembler.dof, neumann_bcs)
     else
       neumann_bcs = NamedTuple()
       neumann_bc_funcs = NamedTuple()
@@ -178,12 +176,12 @@ function Base.show(io::IO, parameters::Parameters)
 end
 
 function create_parameters(
-  assembler, physics, props; 
+  mesh, assembler, physics, props; 
   dirichlet_bcs=DirichletBC[], 
   neumann_bcs=NeumannBC[],
   times=nothing
 )
-  return Parameters(assembler, physics, props, dirichlet_bcs, neumann_bcs, times)
+  return Parameters(mesh, assembler, physics, props, dirichlet_bcs, neumann_bcs, times)
 end
 
 function update_bcs!(p::Parameters)

src/bcs/BoundaryConditions.jl

@@ -31,7 +31,7 @@ end
 $(TYPEDSIGNATURES)
 """
 function BCBookKeeping(
-  dof::DofManager, var_name::Symbol, sset_name::Symbol
+  mesh, dof::DofManager, var_name::Symbol, sset_name::Symbol
 )
   # check if var exists
   # var_index = 0
@@ -66,17 +66,17 @@ function BCBookKeeping(
   fspace = function_space(dof)
 
   # get sset specific fields
-  elements = getproperty(fspace.sideset_elems, sset_name)
-  nodes = getproperty(fspace.sideset_nodes, sset_name)
-  sides = getproperty(fspace.sideset_sides, sset_name)
-  side_nodes = getproperty(fspace.sideset_side_nodes, sset_name)
+  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)
 
   dofs = Vector{Int64}(undef, 0)
   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(fspace.elem_id_maps))
+  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)

src/bcs/DirichletBCs.jl

@@ -44,8 +44,8 @@ $(TYPEDEF)
 $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
-function DirichletBCContainer(dof::DofManager, dbc::DirichletBC)
-  bk = BCBookKeeping(dof, dbc.var_name, dbc.sset_name)
+function DirichletBCContainer(mesh, dof::DofManager, dbc::DirichletBC)
+  bk = BCBookKeeping(mesh, dof, dbc.var_name, dbc.sset_name)
 
   # sort nodes and dofs for dirichlet bc
   dof_perm = _unique_sort_perm(bk.dofs)
@@ -66,10 +66,6 @@ function Base.length(bc::DirichletBCContainer)
   return length(bc.bookkeeping.dofs)
 end
 
-# function bc_set_ids(bc::DirichletBCContainer)
-#   return bc.bookkeeping.nodes
-# end
-
 # need checks on if field types are compatable
 """
 $(TYPEDSIGNATURES)
@@ -143,18 +139,16 @@ function _update_bcs!(bc::DirichletBCContainer, U, backend::KA.Backend)
   return nothing
 end
 
-function create_dirichlet_bcs(dof::DofManager, dirichlet_bcs::Vector{<:DirichletBC})
+function create_dirichlet_bcs(mesh, dof::DofManager, dirichlet_bcs::Vector{<:DirichletBC})
   if length(dirichlet_bcs) == 0
     return NamedTuple(), NamedTuple()
   end
 
   syms = map(x -> Symbol("dirichlet_bc_$x"), 1:length(dirichlet_bcs))
-  # dbcs = NamedTuple{tuple(syms...)}(tuple(dbcs...))
-  # dbcs = DirichletBCContainer(dbcs, size(assembler.dof.H1_vars[1].fspace.coords, 1))
   dirichlet_bc_funcs = NamedTuple{tuple(syms...)}(
     map(x -> x.func, dirichlet_bcs)
   )
-  dirichlet_bcs = DirichletBCContainer.((dof,), dirichlet_bcs)
+  dirichlet_bcs = DirichletBCContainer.((mesh,), (dof,), dirichlet_bcs)
 
   if length(dirichlet_bcs) > 0
     temp_dofs = mapreduce(x -> x.bookkeeping.dofs, vcat, dirichlet_bcs)

src/bcs/NeumannBCs.jl

@@ -96,11 +96,11 @@ end
 
 # TODO below method also currently likely doesn't
 # handle blocks correclty 
-function create_neumann_bcs(dof::DofManager, neumann_bcs::Vector{NeumannBC})
+function create_neumann_bcs(mesh, dof::DofManager, neumann_bcs::Vector{NeumannBC})
   sets = map(x -> x.sset_name, neumann_bcs)
   vars = map(x -> x.var_name, neumann_bcs)
   funcs = map(x -> x.func, neumann_bcs)
-  bks = BCBookKeeping.((dof,), vars, sets)
+  bks = BCBookKeeping.((mesh,), (dof,), vars, sets)
   # bks = _split_bookkeeping_by_block(bks)
   fspace = function_space(dof)
   new_bcs = NeumannBCContainer[]

test/TestAssemblers.jl

@@ -20,7 +20,7 @@ include("poisson/TestPoissonCommon.jl")
     DirichletBC(:u, :sset_3, bc_func),
     DirichletBC(:u, :sset_4, bc_func),
   ]
-  p = create_parameters(asm, physics, props; dirichlet_bcs=dbcs)
+  p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
   Uu = create_unknowns(asm)
   Vu = create_unknowns(asm)
 

test/TestBCs.jl

@@ -14,7 +14,7 @@ function test_dirichlet_bc_container_init()
   u = VectorFunction(fspace, :displ)
   dof = DofManager(u)
   bc_in = DirichletBC(:displ_x, :sset_1, dummy_func_1)
-  bc = FiniteElementContainers.DirichletBCContainer(dof, bc_in)
+  bc = FiniteElementContainers.DirichletBCContainer(mesh, dof, bc_in)
 end
 
 function test_neumann_bc_input()
@@ -30,7 +30,7 @@ function test_neumann_bc_container_init()
   u = VectorFunction(fspace, :displ)
   dof = DofManager(u)
   bc_in = NeumannBC(:displ, :sset_1, dummy_func_2)
-  bc = FiniteElementContainers.NeumannBCContainer(dof, bc_in)
+  bc = FiniteElementContainers.NeumannBCContainer(mesh, dof, bc_in)
 end
 
 @testset ExtendedTestSet "BoundaryConditions" begin

test/mechanics/TestMechanics.jl

@@ -29,7 +29,7 @@ function mechanics_test()
 
   # pre-setup some scratch arrays
   times = TimeStepper(0., 1., 1)
-  p = create_parameters(asm, physics, props; dirichlet_bcs=dbcs, times=times)
+  p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs, times=times)
 
   solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
   integrator = QuasiStaticIntegrator(solver)