allowing dirichlet bcs to be created from either blocks, nodesets, or…

src/FiniteElementContainers.jl

@@ -72,11 +72,11 @@ export update_field_ics!
 export update_ic_values!
 
 # Integrals
-export MatrixIntegral
-export ScalarIntegral
-export VectorIntegral
-export integrate
-export remove_fixed_dofs!
+# export MatrixIntegral
+# export ScalarIntegral
+# export VectorIntegral
+# export integrate
+# export remove_fixed_dofs!
 
 # Integrators
 export AbstractIntegrator

src/assemblers/Constraints.jl

@@ -135,7 +135,7 @@ function _adjust_vector_entries_for_constraints!(b, constraint_storage, ::KA.CPU
     # modify b => (I - G) * b + (Gu - g)
     # but Gu = g, so we don't need that here
     # unless we want to modify this to support weakly
-    # enforced BCs later
+    # enforced BCs later TODO
     for i in 1:length(constraint_storage)
         @inbounds b[i] = (1. - constraint_storage[i]) * b[i]
     end

src/bcs/DirichletBCs.jl

@@ -1,5 +1,7 @@
 abstract type AbstractDirichletBC{F} <: AbstractBC{F} end
 
+const EntityName = Union{Nothing, Symbol}
+
 """
 $(TYPEDEF)
 $(TYPEDSIGNATURES)
@@ -8,7 +10,9 @@ User facing API to define a ```DirichletBC````.
 """
 struct DirichletBC{F} <: AbstractDirichletBC{F}
   func::F
-  sset_name::Symbol
+  block_name::EntityName
+  nset_name::EntityName
+  sset_name::EntityName
   var_name::Symbol
 end
 
@@ -17,17 +21,39 @@ $(TYPEDEF)
 $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
-function DirichletBC(var_name::Symbol, sset_name::Symbol, func::Function)
-  return DirichletBC(func, sset_name, var_name)
+function DirichletBC(
+  var_name::Symbol, func::Function;
+  block_name::EntityName = nothing,
+  nodeset_name::EntityName = nothing,
+  sideset_name::EntityName = nothing
+)
+  return DirichletBC(func, block_name, nodeset_name, sideset_name, var_name)
 end
 
 """
 $(TYPEDEF)
 $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
-function DirichletBC(var_name::String, sset_name::String, func::Function)
-  return DirichletBC(Symbol(var_name), Symbol(sset_name), func)
+function DirichletBC(
+  var_name::String, func::Function;
+  block_name::Union{Nothing, String} = nothing,
+  nodeset_name::Union{Nothing, String} = nothing,
+  sideset_name::Union{Nothing, String} = nothing
+)
+  # return DirichletBC(Symbol(var_name), Symbol(sset_name), func)
+  if block_name !== nothing
+    block_name = Symbol(block_name)
+  end
+
+  if nodeset_name !== nothing
+    nodeset_name = Symbol(nodeset_name)
+  end
+
+  if sideset_name !== nothing
+    sideset_name = Symbol(sideset_name)
+  end
+  return DirichletBC(func, block_name, nodeset_name, sideset_name, Symbol(var_name))
 end
 
 """
@@ -53,7 +79,17 @@ $(TYPEDSIGNATURES)
 $(TYPEDFIELDS)
 """
 function DirichletBCContainer(mesh, dof::DofManager, dbc::DirichletBC)
-  bk = BCBookKeeping(mesh, dof, dbc.var_name, sset_name=dbc.sset_name)
+  if dbc.block_name !== nothing
+    bk = BCBookKeeping(mesh, dof, dbc.var_name, block_name=dbc.block_name)
+  elseif dbc.nset_name !== nothing
+    bk = BCBookKeeping(mesh, dof, dbc.var_name, nset_name=dbc.nset_name)
+  elseif dbc.sset_name !== nothing
+    bk = BCBookKeeping(mesh, dof, dbc.var_name, sset_name=dbc.sset_name)
+  else
+    @assert false
+  end
+
+  # 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)

src/integrals/Integrals.jl

@@ -12,53 +12,86 @@ function Adapt.adapt_storage(to, int::AbstractIntegral)
     )
 end
 
+function Base.show(io::IO, int::AbstractIntegral)
+    println(io, typeof(int).name.name, ":")
+    println(io, "  Assembler type     = $(typeof(int.assembler))")
+    println(io, "  Cache type         = $(typeof(int.cache))")
+    println(io, "  Integrand function = $(int.integrand)")
+end
+
 KA.get_backend(int::AbstractIntegral) = KA.get_backend(int.cache)
 
+function create_field(int::AbstractIntegral)
+    return create_field(int.assembler)
+end
+
+function create_unknowns(int::AbstractIntegral)
+    return create_unknowns(int.assembler)
+end
+
 function integrate end
 
 function (integral::AbstractIntegral)(U, p)
     return integrate(integral, U, p)
 end
 
-struct ScalarIntegral{A, C <: NamedTuple, I} <: AbstractIntegral{A, C, I}
+# abstract type AbstractIntegral{A, C, I} <: AbstractIntegral{A, C, I} end
+# abstract type AbstractBoundaryIntegral{A, C, I} <: AbstractIntegral{A, C, I} end
+
+abstract type AbstractMatrixIntegral{A, C, I} <: AbstractIntegral{A, C, I} end
+abstract type AbstractScalarIntegral{A, C, I} <: AbstractIntegral{A, C, I} end
+abstract type AbstractVectorIntegral{A, C, I} <: AbstractIntegral{A, C, I} end
+
+struct ScalarCellIntegral{A, C <: NamedTuple, I} <: AbstractIntegral{A, C, I}
     assembler::A
     cache::C
     integrand::I
 end
 
-function ScalarIntegral(asm, integrand)
+function ScalarCellIntegral(asm, integrand)
     cache = create_assembler_cache(asm, AssembledScalar())
-    return ScalarIntegral(asm, cache, integrand)
+    return ScalarCellIntegral(asm, cache, integrand)
 end
 
-function gradient(integral::ScalarIntegral)
-    func(physics, interps, x, t, dt, u, u_n, state_old, state_new, props) = ForwardDiff.gradient(
-        z -> integral.integrand(physics, interps, x, t, dt, z, u_n, state_old, state_new, props),
-        u
-    )
+function gradient(integral::ScalarCellIntegral)
+    func(physics, interps, x, t, dt, u, u_n, state_old, state_new, props) = 
+        ForwardDiff.gradient(
+            z -> integral.integrand(physics, interps, x, t, dt, z, u_n, state_old, state_new, props),
+            u
+        )
     cache = create_field(integral.assembler)
-    return VectorIntegral(integral.assembler, cache, func)
+    return VectorCellIntegral(integral.assembler, cache, func)
+end
+
+function hessian(integral::ScalarCellIntegral)
+    func(physics, interps, x, t, dt, u, u_n, state_old, state_new, props) = 
+        ForwardDiff.hessian(
+            z -> integral.integrand(physics, interps, x, t, dt, z, u_n, state_old, state_new, props),
+            u
+        )
+    cache = create_assembler_cache(asm, AssembledMatrix())
+    return MatrixCellIntegral(integral.assembler, cache, func)
 end
 
-function integrate(integral::ScalarIntegral, U, p)
+function integrate(integral::ScalarCellIntegral, U, p)
     cache, dof = integral.cache, integral.assembler.dof
     func = integral.integrand
     assemble_quadrature_quantity!(cache, dof, nothing, func, U, p)
     return mapreduce(sum, sum, values(cache))
 end
 
-struct MatrixIntegral{A, C <: AbstractVector, I} <: AbstractIntegral{A, C, I}
+struct MatrixCellIntegral{A, C <: AbstractVector, I} <: AbstractIntegral{A, C, I}
     assembler::A
     cache::C
     integrand::I
 end
 
-function MatrixIntegral(asm, integrand)
+function MatrixCellIntegral(asm, integrand)
     cache = create_assembler_cache(asm, AssembledMatrix())
-    return MatrixIntegral(asm, cache, integrand)
+    return MatrixCellIntegral(asm, cache, integrand)
 end
 
-function integrate(integral::MatrixIntegral, U, p)
+function integrate(integral::MatrixCellIntegral, U, p)
     assemble_matrix!(
         integral.cache, 
         integral.assembler.matrix_pattern,
@@ -72,7 +105,7 @@ function integrate(integral::MatrixIntegral, U, p)
     )
 end
 
-function remove_fixed_dofs!(integral::MatrixIntegral)
+function remove_fixed_dofs!(integral::MatrixCellIntegral)
     backend = KA.get_backend(integral)
 
     # TODO change API so we don't need to
@@ -87,18 +120,18 @@ function remove_fixed_dofs!(integral::MatrixIntegral)
     return nothing
 end
 
-struct VectorIntegral{A, C <: AbstractField, I} <: AbstractIntegral{A, C, I}
+struct VectorCellIntegral{A, C <: AbstractField, I} <: AbstractIntegral{A, C, I}
     assembler::A
     cache::C
     integrand::I
 end
 
-function VectorIntegral(asm, integrand)
+function VectorCellIntegral(asm, integrand)
     cache = create_assembler_cache(asm, AssembledVector())
-    return VectorIntegral(asm, cache, integrand)
+    return VectorCellIntegral(asm, cache, integrand)
 end
 
-function integrate(integral::VectorIntegral, U, p)
+function integrate(integral::VectorCellIntegral, U, p)
     assemble_vector!(
         integral.cache, 
         integral.assembler.vector_pattern,
@@ -109,7 +142,7 @@ function integrate(integral::VectorIntegral, U, p)
     return integral.cache
 end
 
-function remove_fixed_dofs!(integral::VectorIntegral)
+function remove_fixed_dofs!(integral::VectorCellIntegral)
     backend = KA.get_backend(integral)
     _adjust_vector_entries_for_constraints!(
         integral.cache, integral.assembler.constraint_storage, backend

src/meshes/Exodus.jl

@@ -202,7 +202,7 @@ function PostProcessor(
       append!(element_var_names, names(var))
     # L2QuadratureField case below
     elseif isa(var.fspace.coords, NamedTuple)
-      max_q_points = mapreduce(num_quadrature_points, maximum, values(var.fspace.ref_fes))
+      max_q_points = map(num_quadrature_points, values(var.fspace.ref_fes)) |> maximum
       temp_names = Symbol[]
       for name in names(var)
         for q in 1:max_q_points

src/solvers/DirectLinearSolver.jl

@@ -32,3 +32,30 @@ function solve!(solver::DirectLinearSolver, Uu, p)
   map!((x, y) -> x + y, Uu, Uu, solver.ΔUu)
   return nothing
 end
+# struct DirectLinearSolver{
+#   U <: AbstractArray,
+#   A <: MatrixIntegral,
+#   b <: VectorIntegral,
+#   P
+# }
+#   Δu::U
+#   matrix::A
+#   vector::b
+#   precon::P
+#   timer::TimerOutput
+# end
+
+# function DirectLinearSolver(A::MatrixIntegral, b::VectorIntegral)
+#   precon = I
+#   ΔUu = create_unknowns(A)
+#   return DirectLinearSolver(Δu, A, b, precon, TimerOutput())
+# end
+
+# function solve!(solver::DirectLinearSolver, u, p)
+#   R = solver.matrix(u, p)
+#   K = solver.vector(u, p)
+#   # TODO need bcs as well
+#   solver.Δu .= -K \ R
+#   map!((x, y) -> x + y, u, u, solver.Δu)
+#   return nothing 
+# end

test/TestAssemblers.jl

@@ -26,10 +26,10 @@ function test_sparse_matrix_assembler()
   u = ScalarFunction(V, :u)
   @show asm = SparseMatrixAssembler(u)
   dbcs = DirichletBC[
-    DirichletBC(:u, :sset_1, bc_func),
-    DirichletBC(:u, :sset_2, bc_func),
-    DirichletBC(:u, :sset_3, bc_func),
-    DirichletBC(:u, :sset_4, bc_func),
+    DirichletBC(:u, bc_func; sideset_name = :sset_1),
+    DirichletBC(:u, bc_func; sideset_name = :sset_2),
+    DirichletBC(:u, bc_func; sideset_name = :sset_3),
+    DirichletBC(:u, bc_func; sideset_name = :sset_4),
   ]
   p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
   FiniteElementContainers.initialize!(p)

test/TestBCs.jl

@@ -2,13 +2,49 @@ dummy_func_1(x, t) = 5. * t
 dummy_func_2(x, t) = 5. * SVector{2, Float64}(1., 0.)
 
 function test_dirichlet_bc_input()
-  bc = DirichletBC(:my_var, :my_sset, dummy_func_1)
+  # symbol constructors
+  bc = DirichletBC(:my_var, dummy_func_1; block_name = :my_block)
+  @test bc.block_name == :my_block
+  @test bc.nset_name === nothing
+  @test bc.sset_name === nothing
   @test bc.var_name == :my_var
+  @test typeof(bc.func) == typeof(dummy_func_1)
+
+  bc = DirichletBC(:my_var, dummy_func_1; nodeset_name = :my_nset)
+  @test bc.block_name === nothing
+  @test bc.nset_name == :my_nset
+  @test bc.sset_name === nothing
+  @test bc.var_name == :my_var
+  @test typeof(bc.func) == typeof(dummy_func_1)
+
+
+  bc = DirichletBC(:my_var, dummy_func_1; sideset_name = :my_sset)
+  @test bc.block_name === nothing
+  @test bc.nset_name === nothing
   @test bc.sset_name == :my_sset
+  @test bc.var_name == :my_var
   @test typeof(bc.func) == typeof(dummy_func_1)
-  bc = DirichletBC("my_var", "my_sset", dummy_func_1)
+
+  # string constructors
+  bc = DirichletBC("my_var", dummy_func_1; block_name = "my_block")
+  @test bc.block_name == :my_block
+  @test bc.nset_name === nothing
+  @test bc.sset_name === nothing
+  @test bc.var_name == :my_var
+  @test typeof(bc.func) == typeof(dummy_func_1)
+
+  bc = DirichletBC("my_var", dummy_func_1; nodeset_name = "my_nset")
+  @test bc.block_name === nothing
+  @test bc.nset_name == :my_nset
+  @test bc.sset_name === nothing
   @test bc.var_name == :my_var
+  @test typeof(bc.func) == typeof(dummy_func_1)
+
+  bc = DirichletBC("my_var", dummy_func_1; sideset_name = "my_sset")
+  @test bc.block_name === nothing
+  @test bc.nset_name === nothing
   @test bc.sset_name == :my_sset
+  @test bc.var_name == :my_var
   @test typeof(bc.func) == typeof(dummy_func_1)
 end
 
@@ -17,7 +53,7 @@ function test_dirichlet_bc_container_init()
   fspace = FunctionSpace(mesh, H1Field, Lagrange)
   u = VectorFunction(fspace, :displ)
   dof = DofManager(u)
-  bc_in = DirichletBC(:displ_x, :sset_1, dummy_func_1)
+  bc_in = DirichletBC(:displ_x, dummy_func_1; sideset_name = :sset_1)
   bc = FiniteElementContainers.DirichletBCContainer(mesh, dof, bc_in)
   @show bc
 end

test/TestFunctionSpaces.jl

@@ -88,6 +88,11 @@ end
 #   # # )
 # end
 
+
+function test_bad_interp_type(mesh)
+  @test_throws MethodError FunctionSpace(mesh, H1Field, :SomethingNotSupported)
+end
+
 function test_fspace_h1_field(mesh)
   @show fspace = FunctionSpace(mesh, H1Field, Lagrange)
 end
@@ -109,6 +114,7 @@ function test_function_spaces()
   # coords = H1Field(coords)
 
   mesh = UnstructuredMesh("mechanics/mechanics.g")
+  test_bad_interp_type(mesh)
   test_fspace_h1_field(mesh)
   test_fspace_l2_element_field(mesh)
   test_fspace_l2_quadrature_field(mesh)

test/TestIntegrals.jl

@@ -6,7 +6,7 @@ function test_scalar_integral()
     asm = SparseMatrixAssembler(u)
     f(X, _) = 2. * π^2 * sin(π * X[1]) * sin(π * X[2])
     physics = Poisson(f)
-    integral = FiniteElementContainers.ScalarIntegral(asm, energy)
+    integral = FiniteElementContainers.ScalarCellIntegral(asm, energy)
     grad_integral = FiniteElementContainers.gradient(integral)
 end