diff --git a/src/Decapodes.jl b/src/Decapodes.jl
index 5dbffd89..8b5c28d5 100644
--- a/src/Decapodes.jl
+++ b/src/Decapodes.jl
@@ -24,6 +24,7 @@ export normalize_unicode, DerivOp, append_dot, unicode!, vec_to_dec!,
   CartesianPoint, SpherePoint, r, theta, phi, TangentBasis, θhat, ϕhat,
   average_rewrite, recursive_delete_parents, contract_operators,
   default_dec_matrix_generate, default_dec_generate, 
+  collate,
   @decapode
 
 normalize_unicode(s::String) = Unicode.normalize(s, compose=true, stable=true, chartransform=Unicode.julia_chartransform)
diff --git a/src/collages.jl b/src/collages.jl
new file mode 100644
index 00000000..ba2aa4c8
--- /dev/null
+++ b/src/collages.jl
@@ -0,0 +1,26 @@
+
+"""    function collate(acset, column_names::Vector{Symbol})
+
+Create a collage of two Decapodes that simulates with boundary conditions.
+```
+"""
+function collate(equations, boundaries, uwd, symbols)
+  # TODO: This is assuming only "restriction"-type morphisms.
+
+  f = SummationDecapode{Any, Any, Symbol}()
+  # TODO: Double-check
+  copy_parts!(f, equations, (:Var, :TVar, :Op1, :Op2, :Σ, :Summand, :Sum))
+
+  # TODO: This sets restrictions as Op1s. They are actually Op2s. i.e. Use `bv`.
+  for b in boxes(uwd)
+    ps = incident(uwd, b :box)
+    ev = first(ps)
+    bv = last(ps)
+    var = incident(f, ev, :name)
+    b_var = add_part!(f, :Var, type=f[var, :type], name=f[var, :name])
+    f[var, :name] = Symbol("r_" * string(f[var, :name]))
+    add_part!(f, :Op1, src=b_var, tgt=var, op1=uwd[b, :name])
+  end
+
+  f
+end
diff --git a/test/collages.jl b/test/collages.jl
new file mode 100644
index 00000000..8b2c7a58
--- /dev/null
+++ b/test/collages.jl
@@ -0,0 +1,54 @@
+using Test
+using Decapodes
+using Catlab
+using Catlab.WiringDiagrams
+using Catlab.Programs
+using Catlab.CategoricalAlgebra
+
+
+# TODO: Special names for morphisms that match a method of grabbing boundary
+# simplices.
+
+# TODO: Initial conditions.
+# - This will mean that we need to return both a masked Decapode, and a means of pluggin initial data in for physical quantities, replacing `constuct`.
+# TODO: Temporal boundary conditions.
+# TODO: General boundaries i.e. arbitrary functions of solutions.
+
+# TODO: Add test with empty boundary Decapode.
+
+# Test simple boundary masks.
+DiffusionDynamics = @decapode begin
+  C::Form0
+  ∂ₜ(C) == ∘(d,⋆,d,⋆)(C)
+end
+DiffusionBoundaries = @decapode begin
+  (Cb1, Cb2, Zero)::Form0
+end
+
+DiffusionMorphism = @relation () begin
+  rb1(C, Cb1)
+  rb2(C, Cb2)
+  rb3(Ċ, Zero)
+end
+
+DiffusionCollage = collate(DiffusionMorphism,
+  DiffusionDynamics, DiffusionBoundaries, [
+  (:C, :Cb1),
+  (:C, :Cb2),
+  (:Ċ, :Zero)])
+
+# Note: Since the order does not matter in which rb1 and rb2 are applied, it
+# seems informal to state that one goes before the other.
+# It might be better to provide a semantics for incident edges a la:
+#Diffusion = @decapode begin
+#  C::Form0
+#  #∂ₜ(C) == rb3(∘(d,⋆,d,⋆)(rb1(C)))
+#  #∂ₜ(C) == rb3(∘(d,⋆,d,⋆)(rb2(C)))
+#end
+
+Diffusion = @decapode begin
+  C::Form0
+  ∂ₜ(C) == rb3(∘(d,⋆,d,⋆)(rb2(rb1(C))))
+end
+
+@test_fails is_isomorphic(DiffusionCollage, Diffusion)