diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl
index df42b26..9b17232 100644
--- a/src/FiniteElementContainers.jl
+++ b/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")
diff --git a/src/bcs/BoundaryConditions.jl b/src/bcs/BoundaryConditions.jl
index c4c4703..b667b84 100644
--- a/src/bcs/BoundaryConditions.jl
+++ b/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")
diff --git a/src/bcs/DirichletBCs.jl b/src/bcs/DirichletBCs.jl
index bb46ab3..f9fbcb9 100644
--- a/src/bcs/DirichletBCs.jl
+++ b/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
diff --git a/src/bcs/PeriodicBCs.jl b/src/bcs/PeriodicBCs.jl
new file mode 100644
index 0000000..48b7ab4
--- /dev/null
+++ b/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
diff --git a/src/constraints/Constraints.jl b/src/constraints/Constraints.jl
new file mode 100644
index 0000000..cccaf5d
--- /dev/null
+++ b/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
diff --git a/src/meshes/Meshes.jl b/src/meshes/Meshes.jl
index d27f2d4..58cbada 100644
--- a/src/meshes/Meshes.jl
+++ b/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
diff --git a/src/meshes/StructuredMesh.jl b/src/meshes/StructuredMesh.jl
index b84dfa9..5b0beed 100644
--- a/src/meshes/StructuredMesh.jl
+++ b/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
diff --git a/test/TestBCs.jl b/test/TestBCs.jl
index fe7e998..5bf7d0c 100644
--- a/test/TestBCs.jl
+++ b/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
diff --git a/test/TestMesh.jl b/test/TestMesh.jl
index 73d3187..bd12736 100644
--- a/test/TestMesh.jl
+++ b/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)
 
diff --git a/test/poisson/TestPoissonPBCs.jl b/test/poisson/TestPoissonPBCs.jl
new file mode 100644
index 0000000..3e10847
--- /dev/null
+++ b/test/poisson/TestPoissonPBCs.jl
@@ -0,0 +1,70 @@
+# TODO adapt test to work on GPU
+mesh_file = Base.source_dir() * "/poisson.g"
+
+function test_poisson_pbcs()
+
+    u_analytic(x) = cos(2π*x[1]) + 0.5*cos(4π*x[2]) + 0.25*sin(2π*x[1])*sin(4π*x[2])
+
+    f(X, _) = begin
+        x, y = X
+        (2π)^2 * cos(2π * x) +
+        0.5 * (4π)^2 * cos(4π * y) +
+        0.25 * ((2π)^2 + (4π)^2) * sin(2π * x) * sin(4π * y)
+    end
+    zero_func(_, _) = 0.
+
+    mesh = UnstructuredMesh(mesh_file)
+    V = FunctionSpace(mesh, H1Field, Lagrange) 
+    physics = Poisson(f)
+    props = create_properties(physics)
+    u = ScalarFunction(V, :u)
+    dof = DofManager(u; use_condensed=true)
+    asm = SparseMatrixAssembler(dof)
+    p = create_parameters(mesh, asm, physics, props)
+
+    pbcs = PeriodicBCs(mesh, dof, PeriodicBC[
+        PeriodicBC(:u, :x, :sset_1, :sset_3, zero_func)
+        PeriodicBC(:u, :y, :sset_2, :sset_4, zero_func)
+    ])
+
+    U = create_unknowns(dof)
+    λ = FiniteElementContainers._create_constraint_field(dof, pbcs)
+    sol = vcat(U, λ)
+    assemble_stiffness!(asm, stiffness, U, p)
+
+    K = stiffness(asm)
+    C = FiniteElementContainers._constraint_matrix(dof, pbcs)
+
+    nu = size(K, 1)
+    nc = size(C, 1)
+
+    A = vcat(
+        hcat(K, C'),
+        hcat(C, spzeros(nc, nc))
+    )
+
+    for n in 1:10
+        U = sol[1:nu]
+        λ = sol[nu+1:end]
+        assemble_vector!(asm.residual_storage, asm.dof, residual, U, p)
+        R_u = residual(asm) + C' * λ
+        R_c = C * U
+        R = vcat(R_u, R_c)
+
+        # trying a singel newton step
+        solver = Krylov.MinresSolver(A, -R)
+        Krylov.solve!(solver, A, -R)
+        Δsol = Krylov.solution(solver)
+        sol = sol + Δsol
+
+        # @show norm(R) norm(Δsol)
+        # @show norm(R_u) norm(R_c)
+        # display(λ)
+    end
+
+    u_an = u_analytic.(eachcol(mesh.nodal_coords))
+
+    for n in axes(mesh.nodal_coords, 2)
+        @test isapprox(sol[n], u_an[n], atol=5e-2)
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index cce2499..8067197 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,7 @@ using Aqua
 using CUDA
 using FiniteElementContainers
 using ForwardDiff
+using Krylov
 using LinearAlgebra
 using MPI
 using PartitionedArrays
@@ -12,12 +13,12 @@ using SparseArrays
 using StaticArrays
 using Tensors
 using Test
-# using TestSetExtensions
 
 # put these first so we can use these
 # physics in other tests
 include("mechanics/TestMechanicsCommon.jl")
 include("poisson/TestPoissonCommon.jl")
+include("poisson/TestPoissonPBCs.jl")
 
 include("TestAssemblers.jl")
 include("TestBCs.jl")
@@ -46,6 +47,7 @@ function test_poisson()
       test_poisson(cpu, cond, NewtonSolver, lsolver)
     end
   end
+  test_poisson_pbcs()
 
   if AMDGPU.functional()
     for cond in condensed