diff --git a/docs/Project.toml b/docs/Project.toml
index 6c92ee1..ed35a8e 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,4 @@
 [deps]
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
 FiniteElementContainers = "d08262e4-672f-4e7f-a976-f2cea5767631"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/docs/make.jl b/docs/make.jl
index 7dd68d8..1f3c6e0 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -28,6 +28,13 @@ makedocs(;
     ),
     pages=[
         "Home"                => "index.md",
+        "Tutorial"            => [
+            "1 Poisson Equation"             => "tutorials/1_poisson_equation.md",
+            "2 Advection-Diffusion Equation" => "tutorials/2_advection_diffusion_equation.md",
+            "3 Coupled Problem"              => "tutorials/3_coupled_problem.md",
+            "4 Transient Problem"            => "tutorials/4_transient_problem.md",
+            "5 Solid Mechanics"              => "tutorials/5_solid_mechanics.md"
+        ],
         "Assemblers"          => "assemblers.md",
         "Boundary Conditions" => "boundary_conditions.md",
         "DofManager"          => "dof_manager.md",
diff --git a/docs/src/assemblers.md b/docs/src/assemblers.md
index b099c89..c7991c3 100644
--- a/docs/src/assemblers.md
+++ b/docs/src/assemblers.md
@@ -7,6 +7,13 @@ This section describes the assemblers that are currently available and their abs
 
 All assemblers must possess at minimum a ```DofManager```.
 
+The assemblers possess all the baggage to use the internal method ```sparse!``` in ```SparseArrays.jl```. This method allows for a zero allocation instantiation of a ```SparseMatrixCSC``` type on the CPU. There are also methods available to ease in the conversion of CSC types and other sparse types such as CSR. 
+
+On the GPU, however, this type is first converted to an appropriate COO matrix type on the desired backend. There is unfortunately not a unified sparse matrix API in ```julia``` for GPUs, so we implement this functionality in package extensions. On CUDA, for example, the operational sequence to get a ```CuSparseMatrixCSC``` is to first
+sort the COO ```(row, col, val)``` triplets so they are ordered by row and then column. Then a ```CuSparseMatrixCOO``` type is created and converted to a ```CuSparseMatrixCSC``` type via ```CUDA.jl``` methods. An identical approach is taken for RocM types.
+
+NOTE: This is one of the most actively developed areas of the package. Please use caution with any method beginning with a "_" as these are internal methods that will change without notice.
+
 ## Matrices
 ```@autodocs
 Modules = [FiniteElementContainers]
diff --git a/docs/src/boundary_conditions.md b/docs/src/boundary_conditions.md
index 74301bc..eb3ff24 100644
--- a/docs/src/boundary_conditions.md
+++ b/docs/src/boundary_conditions.md
@@ -8,23 +8,55 @@ and sideset ```sset_1``` with a zero function as follows.
 ```@repl
 using FiniteElementContainers
 bc_func(x, t) = 0.
-bc = DirichletBC(:u, :sset_1, bc_func)
+bc = DirichletBC(:u, bc_func; sideset_name = :sset_1)
 ```
 Internally this is eventually converted in a ```DirichletBCContainer```
 
-### API
+Dirichlet bcs can be setup on element blocks, nodesets, or sidesets. The appropriate keyword argument needs to be supplied with the ```DirichletBC``` constructor.
+
 ```@autodocs
 Modules = [FiniteElementContainers]
 Pages = ["DirichletBCs.jl"]
 Order = [:type, :function]
 ```
 
+# NeumannBC
+We can setup Neumann bcs on a variable ```u``` and sideset ```sset_1``` with a
+simple constant function as follows
+```@repl
+using FiniteElementContainers
+using StaticArrays
+bc_func(x, t) = SVector{1, Float64}(1.)
+bc = NeumannBC(:u, :sset_1, bc_func)
+```
+Note that in comparison to the dirichlet bc example above, the function in this case returns a ```SVector``` of size 1. This will hold for any variable ```u``` that has a single dof. For vector variables, e.g. a traction vector in continuum mechanics, would need something like
+```@repl
+using FiniteElementContainers
+using StaticArrays
+ND = 2
+bc_func(x, t) = SVector{ND, Float64}(1.)
+bc = NeumannBC(:u, :sset_1, bc_func)
+```
+where ```ND``` is the number of dimensions.
+
 ```@autodocs
 Modules = [FiniteElementContainers]
 Pages = ["NeumannBCs.jl"]
 Order = [:type, :function]
 ```
 
+# PeriodicBC
+Periodic boundary conditions are very much a work in progress. There is currently
+some machinary to implement a Lagrange multiplier approach. 
+
+Stay tuned.
+
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["PeriodicBCs.jl"]
+Order = [:type, :function]
+```
+
 ## Boundary Condition Implementation Details
 ```@autodocs
 Modules = [FiniteElementContainers]
diff --git a/docs/src/dof_manager.md b/docs/src/dof_manager.md
index 9a4c17e..7e7b644 100644
--- a/docs/src/dof_manager.md
+++ b/docs/src/dof_manager.md
@@ -18,10 +18,11 @@ mesh = UnstructuredMesh("../../test/poisson/poisson.g")
 V = FunctionSpace(mesh, H1Field, Lagrange)
 u = VectorFunction(V, :u)
 t = ScalarFunction(V, :t)
+f = FiniteElementContainers.GeneralFunction(u, t)
 ```
 Now we can supply these variables to the ```DofManager``` which takes varargs as inputs
 ```@repl dof
-dof = DofManager(u, t)
+dof = DofManager(f)
 ```
 The print methods for this struct show simple metadata about the current dofs for each possible function space.
 
@@ -33,10 +34,10 @@ field = create_unknowns(dof)
 We can create fields of the right size from the ```DofManager``` with the following methods
 
 ```@repl dof
-field = create_field(dof, H1Field)
+field = create_field(dof)
 ```
 
-These methods take the backed of ```dof``` into account to ensure that the fields or unknowns produced are on the same device, e.g. CPU/GPU if ```dof``` is on the CPU/GPU.
+These methods take the backend of ```dof``` into account to ensure that the fields or unknowns produced are on the same device, e.g. CPU/GPU if ```dof``` is on the CPU/GPU.
 
 This struct is created with all dofs initially set as unknown. To modify the unknowns we can do the following
 
@@ -45,4 +46,4 @@ This struct is created with all dofs initially set as unknown. To modify the unk
 Modules = [FiniteElementContainers]
 Pages = ["DofManagers.jl"]
 Order = [:type, :function]
-```
\ No newline at end of file
+```
diff --git a/docs/src/fields.md b/docs/src/fields.md
index d1f53ba..6cf34d4 100644
--- a/docs/src/fields.md
+++ b/docs/src/fields.md
@@ -8,18 +8,18 @@ Fields serve as loose wrappers around ```AbstractArray``` subtypes such that the
 All fields are subtypes of the abstract type ```AbstractField```
 ```@repl
 using FiniteElementContainers
-AbstractField
+FiniteElementContainers.AbstractField
 ```
 
 ## Example - H1Field a.k.a. NodalField
 We can set up a ```H1Field``` in one of two ways. The simplest constructor form can be used as follows
 ```@repl h1field
 using FiniteElementContainers
-field = H1Field(rand(2, 10), (:field_1, :field_2))
+field = H1Field(rand(2, 10))
 ```
 This is stored in a vectorized way as can be seen above
 ```@repl h1field
-field.vals
+field.data
 ```
 Fields can be indexed like regular arrays, e.g.
 ```@repl h1field
@@ -30,11 +30,6 @@ field[1, :]
 ```
 etc.
 
-But they can also be indexed by the symbols provided during construction
-```@repl h1field
-field[:field_1]
-```
-
 ## Abstract type
 The base type for fields is the ```AbstractField``` abstract type. 
 ```@autodocs
diff --git a/docs/src/index.md b/docs/src/index.md
index 6477940..1c6d21c 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -11,16 +11,15 @@ of computational solid mechanics in mind where meshes deform,
 there's path dependence, there's contact between bodies, there are potentially heterogeneous material properties, and other challenges.
 
 If you're primarily interested in writing FEM applications for e.g.
-the Poisson equation or heat equation, there's likely more efficient packages (e.g. ```Gridap``` or ```Ferrite```) out there for this purpose in terms of memory and computational efficiency.
+the Poisson equation or heat equation, there's likely other FEM packages with less mental overhead out there for this purpose.
 
 However, if you need to solve problems with multiple material models, meshes where
-there are mixed elements types, etc. this is likely the only julia package
-at the time of writing this that supports such capabilities. 
+there are mixed elements types, etc. this is the package for you. 
 
 Inspiration for the software design primarily comes from ```fenics``` and ```MOOSE```.
 We've specifically designed the interface to get around all the shortcomings of ```fenics``` 
-(e.g. boundary conditions are a pain, mixed element types a plain, different blocks are pain
+(e.g. boundary conditions are a pain, mixed element types are pain, different blocks are pain
 etc.) 
 
 Our goal is also to ensure all of our methods are next generation hardware capable. This
-means not only supporting things on CPUs but also GPUs (and that doesn't just mean NVIDIA).
\ No newline at end of file
+means not only supporting things on CPUs but also GPUs (and that doesn't just mean NVIDIA). The package is regularly tested against CUDA and RocM aware hardware to ensure all the types, methods, etc. work on CPUs and GPUs. Additionally, we test against Linux, MacOS, and Windows so in theory this should run on any personal computer out of the box. High performance computers are another story.
\ No newline at end of file
diff --git a/docs/src/meshes.md b/docs/src/meshes.md
index c4890bc..da64481 100644
--- a/docs/src/meshes.md
+++ b/docs/src/meshes.md
@@ -3,13 +3,35 @@ CurrentModule = FiniteElementContainers
 ```
 
 # Meshes
-Meshes in ```FiniteElementContainers``` leverage a very abstract interface. No
-single mesh format is directly supported in the main src code but rather different
-mesh types are relegated to package extensions. Currently, only an ```Exodus``` package
-extension is supported but others could be readily supported.
+Meshes in ```FiniteElementContainers``` leverage a very abstract interface. Currently, only an ```Exodus``` interface is directly supported within the main package but others could be readily supported through package extensions which we are planning on.
 
 ```@autodocs
 Modules = [FiniteElementContainers]
 Pages = ["Meshes.jl"]
 Order = [:type, :function]
+```
+
+# Structured Meshes
+Simple structured meshes on rectangles or parallepipeds can be create through ```StructuredMesh``` mesh type.
+
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["StructuredMesh.jl"]
+Order = [:type, :function]
+```
+
+# Unstructured Meshes
+Unstructured meshes (e.g. those read from a file created by a mesher) can be created with the following mesh type
+
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["UnstructuredMesh.jl"]
+Order = [:type, :function]
+```
+
+# Exodus interface API
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["Exodus.jl"]
+Order = [:type, :function]
 ```
\ No newline at end of file
diff --git a/docs/src/tutorials/1_poisson_equation.md b/docs/src/tutorials/1_poisson_equation.md
new file mode 100644
index 0000000..a7ce297
--- /dev/null
+++ b/docs/src/tutorials/1_poisson_equation.md
@@ -0,0 +1,302 @@
+# 1. Poisson Equation
+
+## Strong Form
+Consider the Poisson equation on a domain ``\Omega \subset \mathbb{R}^d`` with boundary ``\partial \Omega``:
+
+``-\nabla \cdot \nabla u = f \quad \text{in } \Omega,``
+ 
+
+with Dirichlet boundary conditions
+
+``u = g \quad \text{on } \partial \Omega_D,``
+
+and Neumann boundary conditions
+
+
+``\nabla u \cdot n = h \quad \text{on } \partial \Omega_N,``
+
+where ``f`` is a source term, and ``n`` is the outward normal vector on the boundary.
+
+## Weak Form
+
+To derive the weak form, multiply the PDE by a test function ``w \in V_0`` which is zero on ``( \partial \Omega_D )`` and integrated over ``\Omega``:
+
+``-\int_\Omega v \, (\nabla \cdot \nabla u) \, d\Omega = \int_\Omega v \, f \, d\Omega``
+
+Applying integration by parts to the left-hand side:
+
+``\int_\Omega \nabla v \cdot \nabla u \, d\Omega - \int_{\partial \Omega} v \, (\nabla u \cdot n) \, d\Gamma = \int_\Omega v \, f \, d\Omega``
+
+Using the boundary conditions:
+
+``\int_\Omega \nabla v \cdot \nabla u \, d\Omega = \int_\Omega v \, f \, d\Omega + \int_{\partial \Omega_N} v \, h \, d\Gamma``
+
+## Finite Element Formulation
+
+Let ``V_h \subset V`` be a finite-dimensional space of trial functions and test functions. The discrete weak form is:
+
+
+``\text{Find } u_h \in V_h \text{ such that } 
+\int_\Omega \nabla v_h \cdot \nabla u_h \, d\Omega = \int_\Omega v_h \, f \, d\Omega + \int_{\partial \Omega_N} v_h \, h \, d\Gamma
+\quad \forall v_h \in V_h.
+``
+
+## Shape function values
+
+Each element ``\Omega_e`` is mapped to a **reference element** ``\hat{\Omega}`` via a geometric mapping:
+
+``
+x(\hat{\xi}) = \sum_{i=1}^{n_\text{nodes}} N_i(\hat{\xi}) \, x_i
+``
+
+where:
+
+- ``N_i(\hat{\xi})`` are the **shape functions** on the reference element  
+- ``x_i`` are the **physical node coordinates**  
+- ``\hat{\xi} \in \hat{\Omega}`` is the reference coordinate  
+
+NOTE: All of this is handled internally. The user does not need to implement this.
+This is just shown for completeness for those not familiar with the nuts and bolts
+of the finite element method.
+
+## Shape function gradients
+
+To compute ``\nabla u_h`` in the physical element, we use the chain rule:
+
+``
+\nabla u_h(x) = \sum_{i=1}^{n_\text{nodes}} u_i \, \nabla N_i(x)
+= \sum_{i=1}^{n_\text{nodes}} u_i \, J^{-T} \, \hat{\nabla} N_i(\hat{\xi})
+``
+
+where:
+
+- ``\hat{\nabla} N_i(\hat{\xi})`` are gradients in the reference element  
+- ``J`` is the **Jacobian matrix** of the mapping:
+
+``
+J = \frac{\partial x}{\partial \hat{\xi}} = \sum_{i=1}^{n_\text{nodes}} x_i \, \hat{\nabla} N_i(\hat{\xi})^T
+``
+
+- ``J^{-T}`` maps reference gradients to physical gradients.  
+
+The determinant of the Jacobian, ``|J|``, scales the quadrature weights when integrating:
+
+NOTE: All of this is handled internally. The user does not need to implement this.
+This is just shown for completeness for those not familiar with the nuts and bolts
+of the finite element method.Shape function values
+
+## Element-wise Assembly with Quadrature
+
+The domain ``\Omega`` is divided into finite elements ``\Omega_e``, and the integrals are computed element-wise:
+
+``
+\int_\Omega \nabla v_h \cdot \nabla u_h \, d\Omega
+= \sum_{e} \int_{\Omega_e} \nabla v_h \cdot \nabla u_h \, d\Omega
+``
+
+Each element integral is approximated using quadrature:
+
+``
+\int_{\Omega_e} \nabla v_h \cdot \nabla u_h \, d\Omega
+\approx \sum_{q=1}^{N_q} w_q \, |J(\hat{\xi}_q)|
+(\nabla v_h(x(\hat{\xi}_q)) \cdot \nabla u_h(x(\hat{\xi}_q)))
+``
+
+where:
+
+- ``x_q`` are the quadrature points in the element ``\Omega_e``  
+- ``w_q`` are the associated quadrature weights  
+- ``N_q`` is the number of quadrature points in the element  
+
+
+## Implementation
+Currently, ```FiniteElementContainers``` expects a set of functions that can be provided to the assemblers that represent
+concepts such as energies, residuals (gradients), and stiffnesses (hessians). Below is an implementation of the above Galerkin finite element formulation in forms suitable to work with the rest of the ```FiniteElementContainers``` machinary. 
+These methods will works on all currently tested backends (e.g. CPU, CUDA, and RocM).
+
+We can take the above equation and represent a discrete residual equation
+
+``
+R_a = \int_\Omega \nabla N^a \cdot \nabla u^a \, d\Omega - \int_\Omega N^a \, f \, d\Omega - \int_{\partial \Omega_N} N^a \, h \, d\Gamma
+``
+than a suitable tangent (stiffness, or whatever you prefer to call it) for a Newton type method is the following
+
+``
+K_{ab} = \frac{\partial R_a}{\partial u_b} = \int_\Omega \nabla N^a \cdot \nabla N^b d\Omega
+``
+
+First we need to create a type for our new physics. We do this by creating a subtype of ```AbstractPhysics``` as follows
+
+```julia
+struct Poisson{F <: Function} <: AbstractPhysics{1, 0, 0}
+  func::F
+end
+```
+
+In the above we create a new type called ```Poisson``` that has a single field, our driving function ``f``. This type is a subtype of ```AbstractPhysics``` which has three
+required generic fields, the number of fields per entity (e.g. node, face, etc.) in the physics (here 1), the number of properties per element, and the number of state variables per quadrature point. For the simple case of a Poisson equation we do not have properties or state variables so these are simply zero. We will explore examples of physical models that require properties and state variables in future tutorials.
+
+Now that we have a type associated with our physics, we can implement the necessary methods for our example. We need a residual method and a stiffness method.
+
+Let's start with the residual method.
+
+```julia
+@inline function FiniteElementContainers.residual(
+  physics::Poisson, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+  interps = map_interpolants(interps, x_el)
+  (; X_q, N, ∇N_X, JxW) = interps
+  ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+  R_q = ∇u_q * ∇N_X' - N' * physics.func(X_q, 0.0)
+  return JxW * R_q[:]
+end
+```
+
+Let's now unpack what probably looks like a fairly complex method at first glance, especially for such a simple equation. 
+- The ```@inline``` macro is required to ensure this runs on GPU backends. Without this you'll likely run into hard to debug GPU errors that may differ from device to device. The first argument to the method is our physics type we defined above. 
+- The second argument is a set of interpolation function values, gradients, and hessians for the reference element of a given block. Note that this method is defined in a such a way that it is element agnostic and operates on a single quadrature point. So the ```interps``` input is a set of shape function values, gradients, hessians, etc. for a given quadrature points. - The next three inputs are the element level coordinates ```x_el```, the current time ```t```, and the current time step ```dt```. 
+- The next two inputs are the current and old element level solution fields ```u_el``` and ```u_el_old```. 
+- These are followed by the quadrature level current and old state variables ```state_old_q``` and ```state_new_q``` and the element level properties ```props_el```. A lot of these fields are completely unnecessary for a problem as simple as the Poisson equation, but they will reveal their usefulness in more compex tutorials. They are required for all physics so we do not need to write a bunch of duplicated assembly methods. Working with varargs on potentially different types on GPUs can be tricky.
+
+Now let's walk through the actual method body. The first two lines
+
+```julia 
+interps = map_interpolants(interps, x_el)
+(; X_q, N, ∇N_X, JxW) = interps
+```
+maps the quadrature level shape function values, and gradients to the current element configuration and unpacks
+the mapped interpolants. ```X_q``` is the quadrature point coordinates in the physical space, ```N``` is the shape function value, ```∇N_X``` is the shape function gradient in the physical space, and ```JxW``` is the determinant of the mapping Jacobian times the Guass quadrature weight.
+
+The next line calls a helper method
+```julia
+∇u_q = interpolate_field_gradients(physics, interps, u_el)
+```
+which uses the shape function gradients and the element level solution field to calculate the quadrature level
+solution field gradient. 
+
+Finally, in the last two lines we compute the residual, integrate it, and return to the calling method.
+```julia
+R_q = ∇u_q * ∇N_X' - N' * physics.func(X_q, 0.0)
+return JxW * R_q[:]
+```
+
+For the stiffness method we can define this analytically as
+```julia
+@inline function FiniteElementContainers.stiffness(
+  physics::Poisson, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+  interps = map_interpolants(interps, x_el)
+  (; X_q, N, ∇N_X, JxW) = interps
+  K_q = ∇N_X * ∇N_X'
+  return JxW * K_q
+end
+```
+or we could use an AD engine if we're feeling lazy.
+
+With these two methods, we have everything we need to assemble the residual and stiffness matrix
+necessary for a standard Newton method. 
+
+## Two-dimensional problem - Poisson on a square
+Let's now solve an actual problem with the methods defined above. Rolling together a simple problem is very straightforward in ```FiniteElementContainers.jl``` and can be achieved in very little code.
+
+Below is an actual example of a homogenous Dirichlet problem on a square. We will use
+the following forcing function 
+
+``
+f(x, y) = 2π^2 sin(πx) sin(πy)
+``
+
+```julia
+using FiniteElementContainers
+
+# setup some helper functions for f and the bcs rhs
+f(X, _) = 2. * π^2 * sin(π * X[1]) * sin(π * X[2])
+bc_func(_, _) = 0.
+
+mesh_file = "my_mesh.exo"
+mesh = UnstructuredMesh(mesh_file)
+V = FunctionSpace(mesh, H1Field, Lagrange) 
+physics = Poisson(f)
+props = create_properties(physics)
+u = ScalarFunction(V, :u)
+asm = SparseMatrixAssembler(u; use_condensed=use_condensed)
+
+# setup bcs
+dbcs = DirichletBC[
+    DirichletBC(:u, bc_func; nodeset_name = :nset_1),
+    DirichletBC(:u, bc_func; nodeset_name = :nset_2),
+    DirichletBC(:u, bc_func; nodeset_name = :nset_3),
+    DirichletBC(:u, bc_func; nodeset_name = :nset_4),
+]
+
+# setup the parameters
+p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
+
+# setup a solver
+solver = NewtonSolver(DirectLinearSolver(asm))
+
+# setup an integrator and let it all evolve one time step
+integrator = QuasiStaticIntegrator(solver)
+evolve!(integrator, p)
+
+# grab our full solution field from our parameters
+U = p.h1_field
+
+# post process results to exodus file
+output_file = "my_output.exo"
+pp = PostProcessor(mesh, output_file, u)
+write_times(pp, 1, 0.0)
+write_field(pp, 1, ("u",), U)
+close(pp)
+```
+
+This has the following visual solution
+
+![](assets/poisson_two_dimensional_dirichlet.png)
+
+# Three dimensional Dirichlet problem
+Now let's look at another problem with same exact code but with a few modifications to boundary conditions, the forcing function, and the mesh file name. We're going to mimic the tutorial from the first tutorial from [MOOSE](https://mooseframework.inl.gov/getting_started/examples_and_tutorials/examples/ex01_inputfile.html)
+
+```julia
+using FiniteElementContainers
+
+# setup some helper functions for f and the bcs rhs
+f(_, _) = 0.0
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+mesh_file = "mug.e"
+mesh = UnstructuredMesh(mesh_file)
+V = FunctionSpace(mesh, H1Field, Lagrange) 
+physics = Poisson(f)
+props = create_properties(physics)
+u = ScalarFunction(V, :u)
+asm = SparseMatrixAssembler(u; use_condensed=use_condensed)
+
+# setup bcs
+dbcs = [
+    DirichletBC(:u, one_func; sideset_name = :bottom)
+    DirichletBC(:u, zero_func; sideset_name = :top)
+]
+# setup the parameters
+p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
+
+# setup a solver
+solver = NewtonSolver(DirectLinearSolver(asm))
+
+# setup an integrator and let it all evolve one time step
+integrator = QuasiStaticIntegrator(solver)
+evolve!(integrator, p)
+
+# grab our full solution field from our parameters
+U = p.h1_field
+
+# post process results to exodus file
+output_file = "my_output.exo"
+pp = PostProcessor(mesh, output_file, u)
+write_times(pp, 1, 0.0)
+write_field(pp, 1, ("u",), U)
+close(pp)
+```
+
+![](assets/poisson_mug.png)
diff --git a/docs/src/tutorials/2_advection_diffusion_equation.md b/docs/src/tutorials/2_advection_diffusion_equation.md
new file mode 100644
index 0000000..aab2671
--- /dev/null
+++ b/docs/src/tutorials/2_advection_diffusion_equation.md
@@ -0,0 +1,113 @@
+# 2. Advection-Diffusion Equation
+
+This example mimics the second [Moose tutorial](https://mooseframework.inl.gov/getting_started/examples_and_tutorials/examples/ex02_kernel.html).
+
+## Strong Form
+Consider the advection equation on a domain ``\Omega \subset \mathbb{R}^d`` with boundary ``\partial \Omega``:
+
+``-\nabla \cdot \nabla u + \mathbf{v}\cdot\nabla u = 0 \quad \text{in } \Omega,``
+ 
+
+with Dirichlet boundary conditions
+
+``u = g \quad \text{on } \partial \Omega_D,``
+
+and Neumann boundary conditions
+
+
+``\nabla u \cdot n = h \quad \text{on } \partial \Omega_N,``
+
+where ``n`` is the outward normal vector on the boundary.
+
+## Weak Form
+
+To derive the weak form, multiply the PDE by a test function ``w \in V_0`` which is zero on ``( \partial \Omega_D )`` and integrated over ``\Omega``:
+
+``-\int_\Omega v \, (\nabla \cdot \nabla u) \, d\Omega + \int_\Omega v \, (\mathbf{v} \cdot \nabla u) \, d\Omega = 0``
+
+Applying integration by parts to the left-hand side:
+
+``\int_\Omega \nabla v \cdot \nabla u \, d\Omega - \int_{\partial \Omega} v \, (\nabla u \cdot n) \, d\Gamma + \int_\Omega v \, (\mathbf{v} \cdot \nabla u) \, d\Omega = 0``
+
+Using the boundary conditions:
+
+``\int_\Omega \nabla v \cdot \nabla u \, d\Omega + \int_\Omega v \, (\mathbf{v} \cdot \nabla u) \, d\Omega = \int_{\partial \Omega_N} v \, h \, d\Gamma``
+
+
+## Implementation
+
+```julia
+struct AdvectionDiffusion{N} <: AbstractPhysics{1, 0, 0}
+    v::SVector{N, Float64}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::AdvectionDiffusion, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u_q = unpack_field(∇u_q, 1)
+    term = dot(∇u_q, physics.v)
+    R_q = ∇N_X * ∇u_q + term * N
+    return JxW * R_q[:]
+end
+  
+@inline function FiniteElementContainers.stiffness(
+    physics::AdvectionDiffusion, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    term = ∇N_X * physics.v
+    K_q = ∇N_X * ∇N_X' + N * term'
+    return JxW * K_q
+end
+```
+
+In the above we gave our new type ```AdvectionDiffusion``` a generic type parameter ```N``` so it can work for 1D, 2D, or 3D problems. It has a single field, the velocity which we could have made a property but since it is constant everywhere, we can just pack it in the type. This physics type has one field and no properties or state variables. The residual and stiffness methods are analogous to the previous example.
+
+# Three dimensional problem
+```julia
+using FiniteElementContainers
+using StaticArrays
+
+# setup some helper functions for f and the bcs rhs
+f(_, _) = 0.0
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+mesh_file = "mug.e"
+mesh = UnstructuredMesh(mesh_file)
+V = FunctionSpace(mesh, H1Field, Lagrange) 
+physics = AdvectionDiffusion(SVector{2, Float64}(0., 0., 1.))
+props = create_properties(physics)
+u = ScalarFunction(V, :u)
+asm = SparseMatrixAssembler(u; use_condensed=use_condensed)
+
+# setup bcs
+dbcs = [
+    DirichletBC(:u, one_func; sideset_name = :bottom)
+    DirichletBC(:u, zero_func; sideset_name = :top)
+]
+# setup the parameters
+p = create_parameters(mesh, asm, physics, props; dirichlet_bcs=dbcs)
+
+# setup a solver
+solver = NewtonSolver(DirectLinearSolver(asm))
+
+# setup an integrator and let it all evolve one time step
+integrator = QuasiStaticIntegrator(solver)
+evolve!(integrator, p)
+
+# grab our full solution field from our parameters
+U = p.h1_field
+
+# post process results to exodus file
+output_file = "my_output.exo"
+pp = PostProcessor(mesh, output_file, u)
+write_times(pp, 1, 0.0)
+write_field(pp, 1, ("u",), U)
+close(pp)
+```
+
+![](assets/advection_diffusion_mug.png)
\ No newline at end of file
diff --git a/docs/src/tutorials/3_coupled_problem.md b/docs/src/tutorials/3_coupled_problem.md
new file mode 100644
index 0000000..452ddce
--- /dev/null
+++ b/docs/src/tutorials/3_coupled_problem.md
@@ -0,0 +1,112 @@
+# 3. Coupled Example
+
+This example mimics the second [Moose tutorial](https://mooseframework.inl.gov/getting_started/examples_and_tutorials/examples/ex03_coupling.html).
+
+## Strong Form
+Consider the advection equation on a domain ``\Omega \subset \mathbb{R}^d`` with boundary ``\partial \Omega``:
+
+``
+\begin{aligned}
+-\nabla \cdot \nabla u + \nabla v \cdot\nabla u &= 0 \quad \text{in } \Omega,\newline
+-\nabla\cdot\nabla v &= 0\quad \text{in } \Omega,
+\end{aligned}
+``
+
+with Dirichlet boundary conditions
+
+``
+\begin{aligned}
+u(x) &= g_1\quad\text{for }x\in\partial\Omega_{d_1} \newline
+v(x) &= g_2\quad\text{for }x\in\partial\Omega_{d_2}
+\end{aligned}
+``
+
+and Neumann boundary conditions
+
+``
+\begin{aligned}
+-\mathbf{n}\cdot\nabla u &= h_1, \newline
+-\mathbf{n}\cdot\nabla v &= h_2
+\end{aligned}
+``
+
+# Weak form
+``
+\begin{aligned}
+\int_\Omega (\nabla w_1\cdot\nabla u + w_1\nabla v\cdot\nabla u)d\Omega &= \int_{\partial\Omega_1} w_1h_1d\Gamma \newline
+\int_\Omega \nabla w_2\cdot\nabla vd\Omega &= \int_{\partial\Omega_2} w_2h_2d\Gamma \newline
+\end{aligned}
+``
+
+# Implementation
+```julia
+struct CoupledPhysics <: AbstractPhysics{2, 0, 0}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::CoupledPhysics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u = unpack_field(∇u_q, 1)
+    ∇v = unpack_field(∇u_q, 2)
+    term = dot(∇u, ∇v)
+    R_u = ∇N_X * ∇u + term * N
+    R_v = ∇N_X * ∇v
+    R = vcat(R_u', R_v')
+    return JxW * R[:]
+end
+
+# begin lazy below
+@inline function FiniteElementContainers.stiffness(
+    physics::CoupledPhysics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    ForwardDiff.jacobian(
+        z -> FiniteElementContainers.residual(
+            physics, interps, x_el, t, dt, z, u_el_old, state_old_q, state_new_q, props_el
+        ), u_el
+    )
+end
+```
+
+# Example problem
+```julia
+one_func(_, _) = 1.0
+two_func(_, _) = 2.0
+zero_func(_, _) = 0.0
+
+mesh = UnstructuredMesh("mug.e")
+V = FunctionSpace(mesh, H1Field, Lagrange)
+physics = CoupledPhysics()
+props = create_properties(physics)
+u = FiniteElementContainers.GeneralFunction(
+    ScalarFunction(V, :u),
+    ScalarFunction(V, :v)
+)
+
+asm = SparseMatrixAssembler(u; use_condensed=true)
+dbcs = [
+    DirichletBC(:u, two_func; sideset_name = :bottom)
+    DirichletBC(:u, zero_func; sideset_name = :top)
+    DirichletBC(:v, one_func; sideset_name = :bottom)
+    DirichletBC(:v, zero_func; sideset_name = :top)
+]
+U = create_field(asm)
+p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
+
+solver = NewtonSolver(DirectLinearSolver(asm))
+integrator = QuasiStaticIntegrator(solver)
+evolve!(integrator, p)
+
+U = p.h1_field
+
+pp = PostProcessor(mesh, "output.e", u)
+write_times(pp, 1, 0.0)
+write_field(pp, 1, ("u", "v"), U)
+close(pp)
+```
+
+Visualized results
+![](assets/coupled_problem_u.png)
+![](assets/coupled_problem_v.png)
diff --git a/docs/src/tutorials/4_transient_problem.md b/docs/src/tutorials/4_transient_problem.md
new file mode 100644
index 0000000..1bcdd25
--- /dev/null
+++ b/docs/src/tutorials/4_transient_problem.md
@@ -0,0 +1,97 @@
+# 4. Transient Problem
+
+This example mimics the second [Moose tutorial](https://mooseframework.inl.gov/getting_started/examples_and_tutorials/examples/ex06_transient.html).
+
+
+## Strong Form
+Consider the Poisson equation on a domain ``\Omega \subset \mathbb{R}^d`` with boundary ``\partial \Omega``:
+
+``\frac{\partial u}{\partial t} -\nabla \cdot \nabla u = 0\quad \text{in } \Omega,``
+ 
+with initial conditions
+``u(x) = u_0(x)``
+
+with Dirichlet boundary conditions
+
+``u = g \quad \text{on } \partial \Omega_D,``
+
+and Neumann boundary conditions
+
+
+``\nabla u \cdot n = h \quad \text{on } \partial \Omega_N,``
+
+## Weak form
+``
+\int_\Omega \nabla v \cdot \nabla u \, d\Omega + \int_\Omega v \frac{\partial u}{\partial t}d\Omega = \int_{\partial \Omega_N} v \, h \, d\Gamma
+``
+where for simplicity we will approximate the time derivate by the following first order approximation
+
+``
+\frac{\partial u}{\partial t} = \frac{u_{n + 1} - u_n}{\Delta t}
+``
+
+## Implementation
+```julia
+struct Transient <: AbstractPhysics{1, 0, 0}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::Transient, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    u_q, ∇u_q = interpolate_field_values_and_gradients(physics, interps, u_el)
+    u_q_old = interpolate_field_values(physics, interps, u_el_old)
+    ∇u = unpack_field(∇u_q, 1)
+    dudt = 20. * (u_q[1] - u_q_old[1]) / dt
+    R = dudt * N + ∇N_X * ∇u
+    return JxW * R[:]
+end
+  
+@inline function FiniteElementContainers.stiffness(
+    physics::Transient, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    K_q = (20. / dt) * N * N' + ∇N_X * ∇N_X'
+    return JxW * K_q
+end
+```
+
+## Example problem
+```julia
+mesh = UnstructuredMesh("cyl-tet.e")
+times = TimeStepper(0., 75., 75)
+V = FunctionSpace(mesh, H1Field, Lagrange)
+physics = Transient()
+props = create_properties(physics)
+u = ScalarFunction(V, :u)
+asm = SparseMatrixAssembler(u; use_condensed=true)
+dbcs = [
+    DirichletBC(:u, zero_func; sideset_name = :bottom)
+    DirichletBC(:u, one_func; sideset_name = :top)
+]
+U = create_field(asm)
+p = create_parameters(
+    mesh, asm, physics, props; 
+    dirichlet_bcs = dbcs,
+    times = times
+)
+solver = NewtonSolver(DirectLinearSolver(asm))
+integrator = QuasiStaticIntegrator(solver)
+pp = PostProcessor(mesh, "output.e", u)
+
+n = 1
+while times.time_current[1] < 75.0
+    evolve!(integrator, p)
+    U = p.h1_field
+    write_times(pp, n, times.time_current[1])
+    write_field(pp, n, ("u",), U)
+    n = n + 1
+end
+
+close(pp)
+```
+
+Visualization
+![](assets/transient.gif)
\ No newline at end of file
diff --git a/docs/src/tutorials/5_solid_mechanics.md b/docs/src/tutorials/5_solid_mechanics.md
new file mode 100644
index 0000000..b2ba98d
--- /dev/null
+++ b/docs/src/tutorials/5_solid_mechanics.md
@@ -0,0 +1,6 @@
+# 5. Solid Mechanics
+For a detailed implementation of solid mechanics
+using ```FiniteElementContainers.jl``` please consult
+the [Cthonios.jl](https://github.com/Cthonios/Cthonios.jl) documentation
+
+TODO add a simple linear elasticity example
diff --git a/docs/src/tutorials/assets/advection_diffusion_mug.png b/docs/src/tutorials/assets/advection_diffusion_mug.png
new file mode 100644
index 0000000..04cdceb
Binary files /dev/null and b/docs/src/tutorials/assets/advection_diffusion_mug.png differ
diff --git a/docs/src/tutorials/assets/coupled_problem_u.png b/docs/src/tutorials/assets/coupled_problem_u.png
new file mode 100644
index 0000000..42b5917
Binary files /dev/null and b/docs/src/tutorials/assets/coupled_problem_u.png differ
diff --git a/docs/src/tutorials/assets/coupled_problem_v.png b/docs/src/tutorials/assets/coupled_problem_v.png
new file mode 100644
index 0000000..2a81a78
Binary files /dev/null and b/docs/src/tutorials/assets/coupled_problem_v.png differ
diff --git a/docs/src/tutorials/assets/poisson_mug.png b/docs/src/tutorials/assets/poisson_mug.png
new file mode 100644
index 0000000..7c8a664
Binary files /dev/null and b/docs/src/tutorials/assets/poisson_mug.png differ
diff --git a/docs/src/tutorials/assets/poisson_two_dimensional_dirichlet.png b/docs/src/tutorials/assets/poisson_two_dimensional_dirichlet.png
new file mode 100644
index 0000000..1b0818e
Binary files /dev/null and b/docs/src/tutorials/assets/poisson_two_dimensional_dirichlet.png differ
diff --git a/docs/src/tutorials/assets/transient.gif b/docs/src/tutorials/assets/transient.gif
new file mode 100644
index 0000000..fe96ed4
Binary files /dev/null and b/docs/src/tutorials/assets/transient.gif differ
diff --git a/examples/moose-tutorials-clone/ex01/script.jl b/examples/moose-tutorials-clone/ex01/script.jl
index 95cc486..6a13cbe 100644
--- a/examples/moose-tutorials-clone/ex01/script.jl
+++ b/examples/moose-tutorials-clone/ex01/script.jl
@@ -40,8 +40,8 @@ function simulate()
     u = ScalarFunction(V, :u)
     asm = SparseMatrixAssembler(u; use_condensed=true)
     dbcs = [
-        DirichletBC(:u, :bottom, one_func)
-        DirichletBC(:u, :top, zero_func)
+        DirichletBC(:u, one_func; sideset_name = :bottom)
+        DirichletBC(:u, zero_func; sideset_name = :top)
     ]
     U = create_field(asm)
     p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
@@ -51,15 +51,15 @@ function simulate()
     # evolve!(integrator, p)
 
     FiniteElementContainers.update_bc_values!(p)
-    residual_int = VectorIntegral(asm, residual)
-    stiffness_int = MatrixIntegral(asm, stiffness)
+    residual_int = FiniteElementContainers.VectorCellIntegral(asm, residual)
+    stiffness_int = FiniteElementContainers.MatrixCellIntegral(asm, stiffness)
 
     K = stiffness_int(U, p)
-    remove_fixed_dofs!(stiffness_int)
+    FiniteElementContainers.remove_fixed_dofs!(stiffness_int)
 
     for n in 1:2
         R = residual_int(U, p)
-        remove_fixed_dofs!(residual_int)
+        FiniteElementContainers.remove_fixed_dofs!(residual_int)
         ΔU = -K \ R.data
         U.data .+= ΔU
         @show n norm(R) norm(ΔU)
diff --git a/examples/moose-tutorials-clone/ex02/script.jl b/examples/moose-tutorials-clone/ex02/script.jl
index 16f63fb..12d85ff 100644
--- a/examples/moose-tutorials-clone/ex02/script.jl
+++ b/examples/moose-tutorials-clone/ex02/script.jl
@@ -3,12 +3,12 @@ using LinearAlgebra
 using StaticArrays
 
 # physics code we actually need to implement
-struct Advection{N} <: AbstractPhysics{1, 0, 0}
+struct AdvectionDiffusion{N} <: AbstractPhysics{1, 0, 0}
     v::SVector{N, Float64}
 end
 
 @inline function FiniteElementContainers.residual(
-    physics::Advection, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+    physics::AdvectionDiffusion, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
 )
     interps = map_interpolants(interps, x_el)
     (; X_q, N, ∇N_X, JxW) = interps
@@ -20,7 +20,7 @@ end
 end
   
 @inline function FiniteElementContainers.stiffness(
-    physics::Advection, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+    physics::AdvectionDiffusion, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
 )
     interps = map_interpolants(interps, x_el)
     (; X_q, N, ∇N_X, JxW) = interps
@@ -37,13 +37,13 @@ function simulate()
     # load a mesh
     mesh = UnstructuredMesh("mug.e")
     V = FunctionSpace(mesh, H1Field, Lagrange)
-    physics = Advection(SVector{3, Float64}(0., 0., 1.))
+    physics = AdvectionDiffusion(SVector{3, Float64}(0., 0., 1.))
     props = create_properties(physics)
     u = ScalarFunction(V, :u)
     asm = SparseMatrixAssembler(u; use_condensed=true)
     dbcs = [
-        DirichletBC(:u, :bottom, one_func)
-        DirichletBC(:u, :top, zero_func)
+        DirichletBC(:u, one_func; sideset_name = :bottom)
+        DirichletBC(:u, zero_func; sideset_name = :top)
     ]
     U = create_field(asm)
     p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
@@ -53,15 +53,15 @@ function simulate()
     # evolve!(integrator, p)
 
     FiniteElementContainers.update_bc_values!(p)
-    residual_int = VectorIntegral(asm, residual)
-    stiffness_int = MatrixIntegral(asm, stiffness)
+    residual_int = FiniteElementContainers.VectorCellIntegral(asm, residual)
+    stiffness_int = FiniteElementContainers.MatrixCellIntegral(asm, stiffness)
 
     K = stiffness_int(U, p)
-    remove_fixed_dofs!(stiffness_int)
+    FiniteElementContainers.remove_fixed_dofs!(stiffness_int)
 
     for n in 1:2
         R = residual_int(U, p)
-        remove_fixed_dofs!(residual_int)
+        FiniteElementContainers.remove_fixed_dofs!(residual_int)
         ΔU = -K \ R.data
         U.data .+= ΔU
         @show n norm(R) norm(ΔU)
diff --git a/examples/moose-tutorials-clone/ex03/script.jl b/examples/moose-tutorials-clone/ex03/script.jl
index 16d584f..00afc45 100644
--- a/examples/moose-tutorials-clone/ex03/script.jl
+++ b/examples/moose-tutorials-clone/ex03/script.jl
@@ -1,3 +1,4 @@
+using AMDGPU
 using FiniteElementContainers
 using ForwardDiff
 using LinearAlgebra
@@ -67,17 +68,20 @@ function simulate()
         ScalarFunction(V, :v)
     )
 
-    asm = SparseMatrixAssembler(u; use_condensed=true)
+    asm = SparseMatrixAssembler(u; use_condensed=false)
     dbcs = [
-        DirichletBC(:u, :bottom, two_func)
-        DirichletBC(:u, :top, zero_func)
-        DirichletBC(:v, :bottom, one_func)
-        DirichletBC(:v, :top, zero_func)
+        DirichletBC(:u, two_func; sideset_name = :bottom)
+        DirichletBC(:u, zero_func; sideset_name = :top)
+        DirichletBC(:v, one_func; sideset_name = :bottom)
+        DirichletBC(:v, zero_func; sideset_name = :top)
     ]
     U = create_field(asm)
     p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
 
-    solver = NewtonSolver(DirectLinearSolver(asm))
+    p = p |> rocm
+    asm = asm |> rocm
+
+    solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
     integrator = QuasiStaticIntegrator(solver)
     evolve!(integrator, p)
 
@@ -96,6 +100,7 @@ function simulate()
     #     @show n norm(R) norm(ΔU)
     # end
 
+    p = p |> cpu
     U = p.h1_field
 
     pp = PostProcessor(mesh, "output.e", u)
diff --git a/examples/moose-tutorials-clone/ex06/script.jl b/examples/moose-tutorials-clone/ex06/script.jl
index 9b57c5d..24a7e80 100644
--- a/examples/moose-tutorials-clone/ex06/script.jl
+++ b/examples/moose-tutorials-clone/ex06/script.jl
@@ -1,5 +1,6 @@
 using FiniteElementContainers
 using LinearAlgebra
+using AMDGPU
 
 # physics code we actually need to implement
 struct Transient <: AbstractPhysics{1, 0, 0}
@@ -41,8 +42,8 @@ function simulate()
     u = ScalarFunction(V, :u)
     asm = SparseMatrixAssembler(u; use_condensed=true)
     dbcs = [
-        DirichletBC(:u, :bottom, zero_func)
-        DirichletBC(:u, :top, one_func)
+        DirichletBC(:u, zero_func; sideset_name = :bottom)
+        DirichletBC(:u, one_func; sideset_name = :top)
     ]
     U = create_field(asm)
     p = create_parameters(
@@ -50,6 +51,11 @@ function simulate()
         dirichlet_bcs = dbcs,
         times = times
     )
+
+    # do it on the GPU whoo!
+    p = p |> rocm
+    asm = asm |> rocm
+
     solver = NewtonSolver(DirectLinearSolver(asm))
     integrator = QuasiStaticIntegrator(solver)
     pp = PostProcessor(mesh, "output.e", u)
@@ -57,7 +63,11 @@ function simulate()
     n = 1
     while times.time_current[1] < 75.0
         evolve!(integrator, p)
-        U = p.h1_field
+
+        # move to cpu for writing
+        # U = p.h1_field
+        p_cpu = p |> cpu
+        U = p_cpu.h1_field
         write_times(pp, n, times.time_current[1])
         write_field(pp, n, ("u",), U)
         n = n + 1
diff --git a/src/DofManagers.jl b/src/DofManagers.jl
index 00fcf6b..375e84c 100644
--- a/src/DofManagers.jl
+++ b/src/DofManagers.jl
@@ -1,8 +1,18 @@
+"""
+$(TYPEDEF)
+$(TYPEDSIGNATURES)
+$(TYPEDFIELDS)
+"""
 abstract type AbstractDofManager{
     IT <: Integer,
     IDs <: AbstractArray{IT, 1}
 } end
 
+"""
+$(TYPEDEF)
+$(TYPEDSIGNATURES)
+$(TYPEDFIELDS)
+"""
 struct DofManager{
     Condensed, # boolean flag for whether or not to seperate data between unknowns and constrained dofs
     # when creating unknowns
@@ -17,6 +27,9 @@ end
 
 # method that initializes dof manager
 # with all dofs unknown
+"""
+$(TYPEDSIGNATURES)
+"""
 function DofManager(var::AbstractFunction; use_condensed::Bool = false)
     dirichlet_dofs = zeros(Int, 0)
     unknown_dofs = 1:size(var.fspace.coords, 2) * length(names(var)) |> collect
@@ -40,14 +53,33 @@ function Adapt.adapt_structure(to, dof::DofManager{C, IT, IDs, Var}) where {C, I
 _field_type(dof::DofManager) = eval(typeof(dof.var.fspace.coords).name.name)
 _is_condensed(dof::DofManager{C, IT, IDs, V}) where {C, IT, IDs, V} = C
 
+"""
+$(TYPEDSIGNATURES)
+Return the total lenth of dofs in the problem.
+
+E.g. for an H1 space this will the number of nodes times
+the number of degrees of freedom per node.
+"""
 Base.length(dof::DofManager) = length(dof.dirichlet_dofs) + length(dof.unknown_dofs)
 
+"""
+$(TYPEDSIGNATURES)
+This returns ```(n_dofs, n_entities)```
+where ```n_dofs``` is the number of dofs per entity (e.g. node)
+and ```n_entities``` is the number of entitities (e.g. node).
+"""
 function Base.size(dof::DofManager)
     l = length(dof)
     nf = length(names(dof.var))
     return (nf, l ÷ nf)
 end
 
+"""
+$(TYPEDSIGNATURES)
+```size(dof, 1)``` returns the number of dofs
+per entity, and ```size(dof, 2)``` returns the 
+number of entities.
+"""
 function Base.size(dof::DofManager, i::Int)
     nf = length(names(dof.var))
     if i == 1
@@ -68,17 +100,42 @@ print(io, "DofManager\n",
 
 KA.get_backend(dof::DofManager) = KA.get_backend(dof.dirichlet_dofs)
 
+"""
+$(TYPEDSIGNATURES)
+Creates a field where ```typeof(field) <: AbstractField```
+based on the variable ```dof``` was created with
+"""
 function create_field(dof::DofManager)
     backend = KA.get_backend(dof)
     field = KA.zeros(backend, Float64, size(dof))
     return _field_type(dof)(field)
 end
 
+"""
+$(TYPEDSIGNATURES)
+Creates a vector of unknown dofs.
+
+This specific method returns a vector equal in length to the
+length of the internally stored list of unknown dofs in ```dof```.
+
+This is used for solution techniques when vector/matrix rows are
+removed where dofs are fixed.
+"""
 function create_unknowns(dof::DofManager{false, IT, IDs, Var}) where {IT, IDs, Var}
     backend = KA.get_backend(dof)
     return KA.zeros(backend, Float64, length(dof.unknown_dofs))
 end
 
+"""
+$(TYPEDSIGNATURES)
+Creates a vector of unknown dofs.
+
+This specific method returns a vector equal in length to the
+length of a field created by ```dof```. E.g. all dofs are unknown.
+
+This is used for solution techniques when vector/matrix rows are not 
+removed where dofs are fixed.
+"""
 function create_unknowns(dof::DofManager{true, IT, IDs, Var}) where {IT, IDs, Var}
     backend = KA.get_backend(dof)
     return KA.zeros(backend, Float64, length(dof))
@@ -116,6 +173,11 @@ function _extract_field_unknowns!(
     return nothing
 end
 
+"""
+$(TYPEDSIGNATURES)
+Updates the entries of ```Uu``` with the unknown dofs
+in the field ```U```.
+"""
 function extract_field_unknowns!(
     Uu::V,
     dof::DofManager{false, IT, IDs, Var},
@@ -132,6 +194,14 @@ function function_space(dof::DofManager)
     return dof.var.fspace
 end
 
+"""
+$(TYPEDSIGNATURES)
+Takes in a list of dof ids associated with dirichlet bcs
+and updates the internals of ```dof``` to reflect these.
+
+NOTE: This clears all existing bcs in ```dof``` and
+starts fresh.
+"""
 function update_dofs!(dof::DofManager, dirichlet_dofs::V) where V <: AbstractArray{<:Integer, 1}
     ND, NI = size(dof)
     Base.resize!(dof.dirichlet_dofs, length(dirichlet_dofs))
@@ -203,6 +273,11 @@ function _update_field_unknowns!(
     return nothing
 end
 
+"""
+$(TYPEDSIGNATURES)
+Updates the unknowns of the field ```U``` based on 
+the values of ```Uu```.
+"""
 function update_field_unknowns!(
     U::AbstractField, 
     dof::DofManager,
@@ -222,6 +297,12 @@ function update_field_unknowns!(
     return nothing
 end
 
+"""
+$(TYPEDSIGNATURES)
+Takes in a field and updates the field
+
+I think this one can be removed/deprecated
+"""
 function update_field_unknowns!(
     U::F, dof::DofManager, Uu::F
 ) where F <: AbstractField
diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl
index bfad648..94f50db 100644
--- a/src/FiniteElementContainers.jl
+++ b/src/FiniteElementContainers.jl
@@ -60,6 +60,7 @@ export num_elements
 # export num_q_points
 
 # Functions
+export GeneralFunction
 export ScalarFunction
 export SymmetricTensorFunction
 export TensorFunction
diff --git a/src/Parameters.jl b/src/Parameters.jl
index c52ed69..8ca9830 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -1,6 +1,16 @@
+"""
+$(TYPEDEF)
+$(TYPEDSIGNATURES)
+$(TYPEDFIELDS)
+"""
 abstract type AbstractParameters end
 
 # TODO need to break up bcs to different field types
+"""
+$(TYPEDEF)
+$(TYPEDSIGNATURES)
+$(TYPEDFIELDS)
+"""
 struct Parameters{
   RT <: Number, # Real type
   RV <: AbstractArray{RT, 1}, # Real vector type
@@ -37,6 +47,9 @@ end
 # 3. figure out how to handle function pointers on the GPU - done
 # 4. add different fspace types
 # 5. convert vectors of dbcs/nbcs into namedtuples - done
+"""
+$(TYPEDSIGNATURES)
+"""
 function Parameters(
   mesh, assembler, physics,
   properties,
@@ -208,6 +221,9 @@ function create_parameters(
   return Parameters(mesh, assembler, physics, props, ics, dirichlet_bcs, neumann_bcs, times)
 end
 
+"""
+$(TYPEDSIGNATURES)
+"""
 function initialize!(p::Parameters)
   update_ic_values!(p.ics, p.h1_coords)
   update_field_ics!(p.h1_field, p.ics)
@@ -230,6 +246,9 @@ function update_bc_values!(p::Parameters)
   return nothing
 end
 
+"""
+$(TYPEDSIGNATURES)
+"""
 function update_dofs!(asm::AbstractAssembler, p::Parameters)
   update_dofs!(asm, p.dirichlet_bcs)
   return nothing
@@ -248,6 +267,9 @@ function _update_for_assembly!(p::Parameters, dof::DofManager, Uu, Vu)
   return nothing
 end
 
+"""
+$(TYPEDSIGNATURES)
+"""
 function update_time!(p::Parameters)
   fill!(p.times.time_current, current_time(p.times) + time_step(p.times))
   return nothing
diff --git a/src/assemblers/SparsityPatterns.jl b/src/assemblers/SparsityPatterns.jl
index 0fc62d3..18d26d2 100644
--- a/src/assemblers/SparsityPatterns.jl
+++ b/src/assemblers/SparsityPatterns.jl
@@ -1,10 +1,10 @@
-"""
-$(TYPEDEF)
-$(TYPEDFIELDS)
-Book-keeping struct for sparse matrices in FEM settings.
-This has all the information to construct a sparse matrix for either
-case where you want to eliminate fixed-dofs or not.
-"""
+# """
+# $(TYPEDEF)
+# $(TYPEDFIELDS)
+# Book-keeping struct for sparse matrices in FEM settings.
+# This has all the information to construct a sparse matrix for either
+# case where you want to eliminate fixed-dofs or not.
+# """
 struct SparseMatrixPattern{
   I <: AbstractArray{Int, 1},
   R <: AbstractArray{Float64, 1}
diff --git a/src/constraints/Constraints.jl b/src/constraints/Constraints.jl
index cccaf5d..a9677d2 100644
--- a/src/constraints/Constraints.jl
+++ b/src/constraints/Constraints.jl
@@ -1,12 +1,12 @@
 abstract type AbstractConstraint end
 abstract type AbstractConstraintContainer end
 
-"""
-Input will be formulated as
+# """
+# Input will be formulated as
 
-c_1 * u_1 + c_2 * u_2 + ... + c_n * u_n = 0
+# c_1 * u_1 + c_2 * u_2 + ... + c_n * u_n = 0
 
-"""
+# """
 struct LinearConstraint{
     RV <: AbstractVector{<:Number},
     RI <: AbstractVector{<:Integer}
diff --git a/src/meshes/StructuredMesh.jl b/src/meshes/StructuredMesh.jl
index 5b0beed..2b8e799 100644
--- a/src/meshes/StructuredMesh.jl
+++ b/src/meshes/StructuredMesh.jl
@@ -1,3 +1,8 @@
+"""
+$(TYPEDEF)
+$(TYPEDSIGNATURES)
+$(TYPEDFIELDS)
+"""
 struct StructuredMesh{
     ND,
     RT <: Number,
diff --git a/src/meshes/UnstructuredMesh.jl b/src/meshes/UnstructuredMesh.jl
index 83ff23c..8506574 100644
--- a/src/meshes/UnstructuredMesh.jl
+++ b/src/meshes/UnstructuredMesh.jl
@@ -37,6 +37,9 @@ function UnstructuredMesh(file_type::AbstractMeshType, file_name::String, create
   return UnstructuredMesh(file, create_edges, create_faces)
 end
 
+"""
+$(TYPEDSIGNATURES)
+"""
 function UnstructuredMesh(file::FileMesh{T}, create_edges::Bool, create_faces::Bool) where T
 
   # read nodal coordinates