periodic-bcs-and-constraints

src/FiniteElementContainers.jl

@@ -20,6 +20,8 @@ export DirichletBC
 export DirichletBCs
 export NeumannBC
 export NeumannBCs
+export PeriodicBC
+export PeriodicBCs
 export update_field_dirichlet_bcs!
 
 # Connectivities
@@ -180,7 +182,9 @@ include("meshes/Meshes.jl")
 include("FunctionSpaces.jl")
 include("Functions.jl")
 include("DofManagers.jl")
+
 include("bcs/BoundaryConditions.jl")
+include("constraints/Constraints.jl")
 include("ics/InitialConditions.jl")
 
 include("formulations/Formulations.jl")

src/bcs/BoundaryConditions.jl

@@ -183,4 +183,5 @@ function update_bc_values!(bcs, funcs, X, t)
 end
 
 include("DirichletBCs.jl")
+include("PeriodicBCs.jl")
 include("NeumannBCs.jl")

src/bcs/DirichletBCs.jl

@@ -221,10 +221,7 @@ end
 function DirichletBCs(mesh, dof, dirichlet_bcs)
 
   if length(dirichlet_bcs) == 0
-    # return NamedTuple(), NamedTuple()
     bc_caches = DirichletBCContainer{Vector{Int}, Vector{Float64}}[]
-    # ti = typeof(identity)
-    # bc_funcs = DirichletBCFunction{ti, ti, ti}[]
     bc_funcs = NamedTuple()
     return DirichletBCs(bc_caches, bc_funcs)
   end
@@ -235,13 +232,9 @@ function DirichletBCs(mesh, dof, dirichlet_bcs)
   )
   dirichlet_bcs = DirichletBCContainer.((mesh,), (dof,), dirichlet_bcs)
 
-  if length(dirichlet_bcs) > 0
-    temp_dofs = mapreduce(x -> x.dofs, vcat, dirichlet_bcs)
-    temp_dofs = unique(sort(temp_dofs))
-    update_dofs!(dof, temp_dofs)
-  end
-
-  # dirichlet_bcs = NamedTuple{tuple(syms...)}(tuple(dirichlet_bcs...))
+  temp_dofs = mapreduce(x -> x.dofs, vcat, dirichlet_bcs)
+  temp_dofs = unique(sort(temp_dofs))
+  update_dofs!(dof, temp_dofs)
 
   return DirichletBCs(dirichlet_bcs, dirichlet_bc_funcs)
 end

src/bcs/PeriodicBCs.jl

@@ -0,0 +1,177 @@
+struct PeriodicBC{F} <: AbstractBC{F}
+    direction::Symbol
+    func::F
+    side_a_sset::Symbol
+    side_b_sset::Symbol
+    var_name::Symbol
+end
+
+function PeriodicBC(
+    var_name::String, direction::String, 
+    side_a_sset::String, side_b_sset::String,
+    func
+)
+    return PeriodicBC(
+        Symbol(direction), func, 
+        Symbol(side_a_sset), Symbol(side_b_sset), 
+        Symbol(var_name)
+    )
+end
+
+function PeriodicBC(
+    var_name::Symbol, direction::Symbol, 
+    side_a_sset::Symbol, side_b_sset::Symbol,
+    func
+)
+    return PeriodicBC(direction, func, side_a_sset, side_b_sset, var_name)
+end
+
+struct PeriodicBCContainer{
+    IV <: AbstractVector{<:Integer},
+    RV <: AbstractVector{<:Number}
+} <: AbstractBCContainer
+    side_a_dofs::IV
+    side_a_nodes::IV
+    side_b_dofs::IV
+    side_b_nodes::IV
+    vals::RV
+end
+
+# TODO add sanity check that tolerance isn't stupid based
+# on mesh sizes
+# TODO also move SVector stuff into barrier func
+function PeriodicBCContainer(
+    mesh, dof::DofManager, pbc::PeriodicBC;
+    tolerance=1.e-6
+)
+    # get direction
+    if pbc.direction == :x
+        dir_id = 1
+    elseif pbc.direction == :y
+        dir_id = 2
+    elseif pbc.direction == :z
+        @assert size(mesh.nodal_coords, 1) == 3
+        dir_id = 3
+    end
+
+    # set up some book keeping
+    side_a_bk = BCBookKeeping(mesh, dof, pbc.var_name, sset_name=pbc.side_a_sset)
+    side_b_bk = BCBookKeeping(mesh, dof, pbc.var_name, sset_name=pbc.side_b_sset)
+
+    # sort the side a sets
+    dof_perm = _unique_sort_perm(side_a_bk.dofs)
+    side_a_dofs = side_a_bk.dofs[dof_perm]
+    side_a_nodes = side_a_bk.nodes[dof_perm]
+    resize!(side_a_bk.dofs, length(side_a_dofs))
+    resize!(side_a_bk.nodes, length(side_a_nodes))
+    copyto!(side_a_bk.dofs, side_a_dofs)
+    copyto!(side_a_bk.nodes, side_a_nodes)
+
+    # now find the nodes in side b associated to those in side a
+    coords = mesh.nodal_coords
+    coords_to_side_a = Dict{Int, Int}()
+    coords_to_side_b = Dict{Int, Int}()
+
+    for node in side_a_nodes
+        temp = round(Int, coords[dir_id, node] / tolerance)
+        coords_to_side_a[temp] = node
+    end
+
+    for node in side_b_bk.nodes
+        temp = round(Int, coords[dir_id, node] / tolerance)
+        coords_to_side_b[temp] = node
+    end
+
+    @assert length(coords_to_side_a) == length(coords_to_side_b) "Side a and side b have different numbers of nodes"
+
+    side_b_nodes = Int[]
+    for node in side_a_nodes
+        temp_coord = round(Int, coords[dir_id, node] / tolerance)
+        side_b_node = coords_to_side_b[temp_coord]
+        push!(side_b_nodes, side_b_node)
+    end
+
+    dofs_all = reshape(1:length(dof), size(dof))
+    dof_id = _dof_index_from_var_name(dof, pbc.var_name)
+    side_b_dofs = Int[]
+    for node in side_b_nodes
+        push!(side_b_dofs, dofs_all[dof_id, node])
+    end
+
+    vals = zeros(length(side_a_dofs))
+
+    return PeriodicBCContainer(
+        side_a_dofs, side_a_nodes,
+        side_b_dofs, side_b_nodes,
+        vals
+    )
+end
+
+function Adapt.adapt_structure(to, bc::PeriodicBCContainer)
+    return PeriodicBCContainer(
+        adapt(to, bc.side_a_dofs),
+        adapt(to, bc.side_a_nodes),
+        adapt(to, bc.side_b_dofs),
+        adapt(to, bc.side_b_nodes),
+        adapt(to, bc.vals)
+    )
+end
+
+struct PeriodicBCFunction{F} <: AbstractBCFunction{F}
+    func::F
+end
+
+struct PeriodicBCs{
+    IV      <: AbstractVector{<:Integer},
+    RV      <: AbstractVector{<:Number},
+    BCFuncs <: NamedTuple
+}
+    bc_caches::Vector{PeriodicBCContainer{IV, RV}}
+    bc_funcs::BCFuncs
+end
+
+function PeriodicBCs(mesh, dof, periodic_bcs)
+    if length(periodic_bcs) == 0
+        bc_caches = PeriodicBCContainer{Vector{Int}, Vector{Float64}}[]
+        bc_funcs = NamedTuple()
+        return PeriodicBCs(bc_caches, bc_funcs)
+    end
+
+    syms = map(x -> Symbol("periodic_bc_$x"), 1:length(periodic_bcs))
+    periodic_bc_funcs = NamedTuple{tuple(syms...)}(
+        map(x -> PeriodicBCFunction(x.func), periodic_bcs)
+    )
+    periodic_bcs = PeriodicBCContainer.((mesh,), (dof,), periodic_bcs)
+
+    # TODO add hooks for dof removal if that's the way we
+    # want to implement it downstream, but make it an optional
+    # kwarg
+
+    return PeriodicBCs(periodic_bcs, periodic_bc_funcs)
+end
+
+function _constraint_matrix(dof::DofManager, pbcs::PeriodicBCs)
+    Is, Js, Vs = Int[], Int[], Float64[]
+    n = 1
+    for bc in pbcs.bc_caches
+        for (dof_a, dof_b) in zip(bc.side_a_dofs, bc.side_b_dofs)
+            # side a dof contribution
+            push!(Is, n)
+            push!(Js, dof_a)
+            push!(Vs, -1.)
+
+            # side b dof contribution
+            push!(Is, n)
+            push!(Js, dof_b)
+            push!(Vs, 1.)
+            n = n + 1
+        end
+    end
+
+    return sparse(Is, Js, Vs)
+end
+
+function _create_constraint_field(dof::DofManager, pbcs::PeriodicBCs)
+    n_constraints = mapreduce(x -> length(x.side_a_dofs), +, pbcs.bc_caches)
+    return zeros(n_constraints)
+end

src/constraints/Constraints.jl

@@ -0,0 +1,29 @@
+abstract type AbstractConstraint end
+abstract type AbstractConstraintContainer end
+
+"""
+Input will be formulated as
+
+c_1 * u_1 + c_2 * u_2 + ... + c_n * u_n = 0
+
+"""
+struct LinearConstraint{
+    RV <: AbstractVector{<:Number},
+    RI <: AbstractVector{<:Integer}
+} <: AbstractConstraint
+    coefficients::RV
+    dofs::RI
+
+    function LinearConstraint(coefficients, dofs)
+        @assert length(coefficients) == length(dofs)
+        new{typeof(coefficients), typeof(dofs)}(coefficients, dofs)
+    end
+end
+
+function evaluate(cons::LinearConstraint, u)
+    val = zero(eltype(u))
+    for n in axes(cons.dofs, 1)
+        val += cons.coefficients[n] * u[cons.dofs[n]]
+    end
+    return val
+end

src/meshes/Meshes.jl

@@ -24,6 +24,31 @@ $(TYPEDEF)
 """
 abstract type AbstractMesh end
 
+function Base.show(io::IO, mesh::AbstractMesh)
+  println(io, typeof(mesh).name.name, ":")
+  println(io, "  Number of dimensions = $(size(mesh.nodal_coords, 1))")
+  println(io, "  Number of nodes      = $(size(mesh.nodal_coords, 2))")
+
+  println(io, "  Element Blocks:")
+  for (conn, name, type) in zip(mesh.element_conns, mesh.element_block_names, mesh.element_types)
+    println(io, "    $name:")
+    println(io, "      Element type       = $type")
+    println(io, "      Number of elements = $(size(conn, 2))")
+  end
+
+  println(io, "  Node sets:")
+  for (name, nodes) in mesh.nodeset_nodes
+    println(io, "    $name:")
+    println(io, "      Number of nodes = $(length(nodes))")
+  end
+
+  println(io, "  Side sets:")
+  for (name, sides) in mesh.sideset_sides
+    println(io, "    $name:")
+    println(io, "      Number of elements = $(length(sides))")
+  end
+end
+
 """
 $(TYPEDSIGNATURES)
 Returns file name for an mesh type

src/meshes/StructuredMesh.jl

@@ -71,15 +71,6 @@ function StructuredMesh(element_type, mins, maxs, counts)
     )
 end
 
-function Base.show(io::IO, mesh::StructuredMesh)
-    println(io, "StructuredMesh:")
-    println(io, "  Number of dimensions = $(size(mesh.nodal_coords, 1))")
-    println(io, "  Number of nodes      = $(size(mesh.nodal_coords, 2))")
-    println(io, "  Number of elements   = $(size(mesh.element_conns[:block_1], 2))")
-    println(io, "  Element type         = $(mesh.element_types[1])")
-
-end
-
 function _hex8_structured_mesh_data(mins, maxs, counts)
     N_x, N_y, N_z = counts
     E_x = N_x - 1

test/TestBCs.jl

@@ -43,9 +43,25 @@ function test_neumann_bc_container_init()
   @show bc
 end
 
+function test_periodic_bc_input()
+  bc = PeriodicBC(:my_var, :x, :my_sset_1, :my_sset_2, dummy_func_1)
+  @test bc.var_name == :my_var
+  @test bc.direction == :x
+  @test bc.side_a_sset == :my_sset_1
+  @test bc.side_b_sset == :my_sset_2
+  @test typeof(bc.func) == typeof(dummy_func_1)
+  bc = PeriodicBC("my_var", "x", "my_sset_1", "my_sset_2", dummy_func_1)
+  @test bc.var_name == :my_var
+  @test bc.direction == :x
+  @test bc.side_a_sset == :my_sset_1
+  @test bc.side_b_sset == :my_sset_2
+  @test typeof(bc.func) == typeof(dummy_func_1)
+end
+
 @testset "BoundaryConditions" begin
   test_dirichlet_bc_input()
   test_dirichlet_bc_container_init()
   test_neumann_bc_input()
   # test_neumann_bc_container_init()
+  test_periodic_bc_input()
 end

test/TestMesh.jl

@@ -68,7 +68,8 @@ function test_structured_mesh()
 end
 
 function test_write_mesh()
-  mesh = StructuredMesh("hex", (0., 0., 0.), (1., 1., 1.), (3, 3, 3))
+  # show below is to cover Base.show in tests
+  @show mesh = StructuredMesh("hex", (0., 0., 0.), (1., 1., 1.), (3, 3, 3))
   FiniteElementContainers.write_to_file(mesh, "shex.exo")
   rm("shex.exo", force = true)