diff --git a/Project.toml b/Project.toml
index 0bdfc13..0e38a00 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
-version = "0.13.2"
+version = "0.13.3"
 authors = ["Craig M. Hamel <cmhamel32@gmail.com> and contributors"]
 
 [deps]
diff --git a/docs/make.jl b/docs/make.jl
index 1f3c6e0..4b46561 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -35,9 +35,11 @@ makedocs(;
             "4 Transient Problem"            => "tutorials/4_transient_problem.md",
             "5 Solid Mechanics"              => "tutorials/5_solid_mechanics.md"
         ],
+        "AppTools"            => "app_tools.md",
         "Assemblers"          => "assemblers.md",
         "Boundary Conditions" => "boundary_conditions.md",
         "DofManager"          => "dof_manager.md",
+        "Expressions"         => "expressions.md",
         "Fields"              => "fields.md",
         "Formulations"        => "formulations.md",
         "Function spaces"     => "function_spaces.md",
diff --git a/docs/src/app_tools.md b/docs/src/app_tools.md
new file mode 100644
index 0000000..71bd9fc
--- /dev/null
+++ b/docs/src/app_tools.md
@@ -0,0 +1,6 @@
+# AppTools
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["AppTools.jl"]
+Order = [:module, :type, :function]
+```
\ No newline at end of file
diff --git a/docs/src/assemblers.md b/docs/src/assemblers.md
index c7991c3..b1abde7 100644
--- a/docs/src/assemblers.md
+++ b/docs/src/assemblers.md
@@ -14,6 +14,13 @@ sort the COO ```(row, col, val)``` triplets so they are ordered by row and then
 
 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.
 
+## Matrix Diagonal
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["Diagonal.jl"]
+Order = [:function]
+```
+
 ## Matrices
 ```@autodocs
 Modules = [FiniteElementContainers]
@@ -42,6 +49,13 @@ Pages = ["QuadratureQuantity.jl"]
 Order = [:function]
 ```
 
+## Sources
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["Source.jl"]
+Order = [:function]
+```
+
 ## Vector
 ```@autodocs
 Modules = [FiniteElementContainers]
@@ -49,6 +63,13 @@ Pages = ["Vector.jl"]
 Order = [:function]
 ```
 
+## Weakly enforced BCs
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["WeaklyEnforcedBCs.jl"]
+Order = [:function]
+```
+
 ## Abstract Interface
 ```@autodocs
 Modules = [FiniteElementContainers]
@@ -56,6 +77,13 @@ Pages = ["Assemblers.jl"]
 Order = [:type, :function]
 ```
 
+## Assembly methods that are safe with Enzyme.jl
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["Enzyme.jl"]
+Order = [:type, :function]
+```
+
 ## SparseMatrixAssembler
 ```@autodocs
 Modules = [FiniteElementContainers]
diff --git a/docs/src/boundary_conditions.md b/docs/src/boundary_conditions.md
index eb3ff24..096e791 100644
--- a/docs/src/boundary_conditions.md
+++ b/docs/src/boundary_conditions.md
@@ -8,7 +8,7 @@ and sideset ```sset_1``` with a zero function as follows.
 ```@repl
 using FiniteElementContainers
 bc_func(x, t) = 0.
-bc = DirichletBC(:u, bc_func; sideset_name = :sset_1)
+bc = DirichletBC("u", bc_func; sideset_name = "sset_1")
 ```
 Internally this is eventually converted in a ```DirichletBCContainer```
 
@@ -27,7 +27,7 @@ simple constant function as follows
 using FiniteElementContainers
 using StaticArrays
 bc_func(x, t) = SVector{1, Float64}(1.)
-bc = NeumannBC(:u, :sset_1, bc_func)
+bc = NeumannBC("u", bc_func, "sset_1")
 ```
 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
@@ -35,7 +35,7 @@ using FiniteElementContainers
 using StaticArrays
 ND = 2
 bc_func(x, t) = SVector{ND, Float64}(1.)
-bc = NeumannBC(:u, :sset_1, bc_func)
+bc = NeumannBC("u", bc_func, "sset_1")
 ```
 where ```ND``` is the number of dimensions.
 
@@ -57,6 +57,20 @@ Pages = ["PeriodicBCs.jl"]
 Order = [:type, :function]
 ```
 
+# RobinBC
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["RobinBCs.jl"]
+Order = [:type, :function]
+```
+
+# Source
+```@autodocs
+Modules = [FiniteElementContainers]
+Pages = ["Sources.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 7e7b644..6cfc312 100644
--- a/docs/src/dof_manager.md
+++ b/docs/src/dof_manager.md
@@ -13,11 +13,11 @@ The ```DofManager``` is a struct that keeps track of which dofs are unknown or c
 A ```DofManager``` can be created as follows. First we must create functions for our variables of interest from their associated function spaces.
 
 ```@repl dof
-using Exodus, FiniteElementContainers
+using FiniteElementContainers
 mesh = UnstructuredMesh("../../test/poisson/poisson.g")
 V = FunctionSpace(mesh, H1Field, Lagrange)
-u = VectorFunction(V, :u)
-t = ScalarFunction(V, :t)
+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
diff --git a/docs/src/expressions.md b/docs/src/expressions.md
new file mode 100644
index 0000000..84b25c3
--- /dev/null
+++ b/docs/src/expressions.md
@@ -0,0 +1,33 @@
+# Expressions
+This package supports type-stable expression functions that can be parsed from strings at runtime. This enables arbitrary boundary condition support in compiled binaries. This feature heavily leverages the package [DynamicExpressions.jl](https://github.com/SymbolicML/DynamicExpressions.jl)
+
+## ScalarExpressionFunction
+```@repl
+import FiniteElementContainers.Expressions: ScalarExpressionFunction
+using StaticArrays
+expr = "5.0 * exp(-2 * t) * cos(2 * pi * t)"
+var_names = ["x", "y", "z", "t"]
+func = ScalarExpressionFunction{Float64}(expr, var_names)
+X = SVector{3, Float64}(1., 2., 3.)
+t = 5.0
+val = func(X, t)
+```
+
+## VectorExpressionFunction
+```@repl
+import FiniteElementContainers.Expressions: VectorExpressionFunction
+using StaticArrays
+exprs = ["0", "5.0 * exp(-2 * t) * cos(2 * pi * t)", "0"]
+var_names = ["x", "y", "z", "t"]
+func = VectorExpressionFunction{3, Float64}(exprs, var_names)
+X = SVector{3, Float64}(1., 2., 3.)
+t = 5.0
+val = func(X, t)
+```
+
+# API
+```@autodocs
+Modules = [FiniteElementContainers.Expressions]
+Pages = ["Expressions.jl"]
+Order = [:type, :function]
+```
\ No newline at end of file
diff --git a/docs/src/fields.md b/docs/src/fields.md
index 36f8ef8..49d0d79 100644
--- a/docs/src/fields.md
+++ b/docs/src/fields.md
@@ -30,57 +30,10 @@ field[1, :]
 ```
 etc.
 
-## Abstract type
+## Fields
 The base type for fields is the ```AbstractField``` abstract type. 
 ```@autodocs
 Modules = [FiniteElementContainers]
-Pages = ["fields/Fields.jl"]
-Order = [:type]
-```
-Any new field added to ```FiniteElementContainers``` should be a subtype of this type.
-
-## Methods for ```AbstractField```
-```@autodocs
-Modules = [FiniteElementContainers]
-Pages = ["fields/Fields.jl"]
-Order = [:function]
-```
-
-## Implementations
-The existing direct subtypes of ```AbstractField``` are the following
-
-### Connectivity
-The connectivity type is a simple alias for ```L2ElementField``` defined below
-```@autodocs
-Modules = [FiniteElementContainers]
-Pages = ["fields/Connectivity.jl"]
-Order = [:type, :function]
-```
-
-### H1 field
-```@autodocs
-Modules = [FiniteElementContainers]
-Pages = ["fields/H1Field.jl"]
-Order = [:type, :function]
-```
-
-### Hcurl field
-```@autodocs
-Modules = [FiniteElementContainers]
-Pages = ["fields/HcurlField.jl"]
+Pages = ["Fields.jl"]
 Order = [:type, :function]
-```
-
-### Hdiv field
-```@autodocs
-Modules = [FiniteElementContainers]
-Pages = ["fields/HdivField.jl"]
-Order = [:type, :function]
-```
-
-### L2 field
-```@autodocs
-Modules = [FiniteElementContainers]
-Pages = ["fields/L2Field.jl"]
-Order = [:type, :function]
-```
+```
\ No newline at end of file
diff --git a/docs/src/formulations.md b/docs/src/formulations.md
index 44002b7..f837501 100644
--- a/docs/src/formulations.md
+++ b/docs/src/formulations.md
@@ -4,7 +4,7 @@ CurrentModule = FiniteElementContainers
 
 # Formulations
 ```@docs
-AbstractMechanicsFormulation
+AbstractElementFormulation
 ```
 
 # Implementations
@@ -12,10 +12,6 @@ AbstractMechanicsFormulation
 PlaneStrain
 ```
 
-```@docs
-ScalarFormulation
-```
-
 ```@docs
 ThreeDimensional
 ```
@@ -29,4 +25,8 @@ extract_stiffness
 extract_stress
 modify_field_gradients
 num_dimensions
+scatter_with_gradients!
+scatter_with_gradients_and_gradients!
+scatter_with_values!
+scatter_with_values_and_values!
 ```
diff --git a/examples/app/input-file.toml b/examples/app/input-file.toml
index 74ce289..f908d53 100644
--- a/examples/app/input-file.toml
+++ b/examples/app/input-file.toml
@@ -1,31 +1,34 @@
 [device]
 backend = "cpu"
 
-[functions.zero_ic]
-type = "scalar expression"
-expression = "0.0"
-variables = ["x", "y"]
-
 [functions.zero_bc]
 type = "scalar expression"
 expression = "0.0"
 variables = ["x", "y", "t"]
 
-[functions.zero_neumann]
+[functions.source_func]
 type = "vector expression"
-expressions = ["1.0"]
+expressions = ["2 * pi^2 * sin(2 * pi * x) * sin(2 * pi * y)"]
 variables = ["x", "y", "t"]
 
 [mesh]
 file_path = "poisson.g"
 file_type = "exodus"
 
-[[initial_conditions]]
-function = "zero_ic"
-blocks = ["block_1"]
+[[boundary_conditions.dirichlet]]
+function = "zero_bc"
+side_sets = ["sset_1", "sset_2"]
 variables = ["u"]
 
 [[boundary_conditions.dirichlet]]
 function = "zero_bc"
-side_sets = ["sset_1", "sset_2", "sset_3", "sset_4"]
+side_sets = ["sset_3", "sset_4"]
 variables = ["u"]
+
+[[boundary_conditions.source]]
+function = "source_func"
+blocks = ["block_1"]
+variables = ["u"]
+
+# [linear_solver.krylov]
+
diff --git a/examples/app/src/MyApp.jl b/examples/app/src/MyApp.jl
index db38e96..bea344b 100644
--- a/examples/app/src/MyApp.jl
+++ b/examples/app/src/MyApp.jl
@@ -9,7 +9,7 @@ using TimerOutputs
 # include files
 include("Physics.jl")
 
-f(X, _) = 2. * π^2 * sin(2π * X[1]) * sin(2π * X[2])
+# f(X, _) = 2. * π^2 * sin(2π * X[1]) * sin(2π * X[2])
 
 const N = 1
 
@@ -38,10 +38,12 @@ function app_main(ARGS::Vector{String})
     asm = SparseMatrixAssembler{SPT, false, false}(dof)
 
     # setup physics and properties
-    physics = Poisson{typeof(f)}[]
+    # physics = Poisson{typeof(f)}[]
+    physics = Laplace[]
     props = Vector{Float64}[]
     for _ in 1:length(sim.mesh.element_conns)
-        temp_physics = Poisson(f)
+        # temp_physics = Poisson(f)
+        temp_physics = Laplace()
         push!(physics, temp_physics)
         push!(props, create_properties(temp_physics))
     end
diff --git a/examples/app/src/Physics.jl b/examples/app/src/Physics.jl
index b3e29fc..52cf8ca 100644
--- a/examples/app/src/Physics.jl
+++ b/examples/app/src/Physics.jl
@@ -1,3 +1,25 @@
+struct Laplace <: AbstractPhysics{1, 0, 0}
+end
+
+@inline function FiniteElementContainers.residual(
+  physics::Laplace, 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'
+  return JxW * R_q[:]
+end
+
+@inline function FiniteElementContainers.stiffness(
+  physics::Laplace, 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
+
 struct Poisson{F <: Function} <: AbstractPhysics{1, 0, 0}
     func::F
 end
diff --git a/src/AppTools.jl b/src/AppTools.jl
index 7020d4c..02e56cb 100644
--- a/src/AppTools.jl
+++ b/src/AppTools.jl
@@ -1,8 +1,3 @@
-"""
-All of these tools are meant to be type safe
-from the perspective of working with JuliaC
-trim mode for small binaries
-"""
 module AppTools
 
 export App
@@ -406,7 +401,7 @@ struct BCSettings{N, T <: Number}
                 vars = temp["variables"]::Vector{String}
                 for side_set in sidesets
                     for var in vars
-                        push!(nbcs, RobinBC(var, func, side_set))
+                        push!(rbcs, RobinBC(var, func, side_set))
                     end
                 end
             end
@@ -419,11 +414,11 @@ struct BCSettings{N, T <: Number}
                 blocks = temp["blocks"]::Vector{String}
                 func = temp["function"]::String
                 func = functions.vector_expr_funcs[func]
-                sidesets = temp["side_sets"]::Vector{String}
+                sidesets = temp["blocks"]::Vector{String}
                 vars = temp["variables"]::Vector{String}
                 for block in blocks
                     for var in vars
-                        push!(nbcs, Source(var, func, block))
+                        push!(srcs, Source(var, func, block))
                     end
                 end
             end
@@ -444,17 +439,19 @@ struct ICSettings{T <: Number}
 
     function ICSettings{T}(log_file, data, functions::FunctionSettings{N, T}) where {N, T}
         print_banner(log_file, "Initial conditions")
-        ic_settings = data["initial_conditions"]::Vector{Any}
         ics = InitialCondition{ScalarExpressionFunction{Float64}}[]
-        for ic in ic_settings
-            temp = ic::Dict{String, Any}
-            blocks = temp["blocks"]::Vector{String}
-            func = temp["function"]::String
-            func = functions.scalar_expr_funcs[func]
-            vars = temp["variables"]::Vector{String}
-            for block in blocks
-                for var in vars
-                    push!(ics, InitialCondition(var, func, block))
+        if haskey(data, "initial_conditions")
+            ic_settings = data["initial_conditions"]::Vector{Any}
+            for ic in ic_settings
+                temp = ic::Dict{String, Any}
+                blocks = temp["blocks"]::Vector{String}
+                func = temp["function"]::String
+                func = functions.scalar_expr_funcs[func]
+                vars = temp["variables"]::Vector{String}
+                for block in blocks
+                    for var in vars
+                        push!(ics, InitialCondition(var, func, block))
+                    end
                 end
             end
         end
@@ -485,6 +482,7 @@ struct InputSettings{N, T <: Number}
     functions::FunctionSettings{N, T}
     ics::ICSettings{T}
     mesh::MeshSettings
+    raw_input::Dict{String, Any}
 end
 
 function InputSettings{N}(cli_args::CLIArgParser, log_file::LogFile, ::Type{T} = Float64) where {N, T <: Number}
@@ -503,7 +501,7 @@ function InputSettings{N}(cli_args::CLIArgParser, log_file::LogFile, ::Type{T} =
     bcs = BCSettings{N, T}(log_file, data, functions)
     ics = ICSettings{T}(log_file, data, functions)
     mesh = MeshSettings(log_file, data)
-    return InputSettings{N, T}(bcs, functions, ics, mesh)
+    return InputSettings{N, T}(bcs, functions, ics, mesh, data)
 end
 
 # function _parse_function_space_settings(log_file, data)
diff --git a/src/Expressions.jl b/src/Expressions.jl
index 644e09d..194ba3f 100644
--- a/src/Expressions.jl
+++ b/src/Expressions.jl
@@ -1,6 +1,7 @@
 module Expressions
 
 import DynamicExpressions.NodeModule: DEFAULT_MAX_DEGREE
+using DocStringExtensions
 using DynamicExpressions
 using StaticArrays
 
@@ -307,6 +308,8 @@ function _nud(p::Parser, t::Token, ::Type{T}) where T <: Number
             if string(t.name) in p.var_names
                 name_id = findfirst(x -> x == string(t.name), p.var_names)
                 return Node{T}(; feature = name_id)
+            elseif string(t.name) == "pi"
+                return Node{T}(; val = π)
             # elseif string(t.name) in p.parameter_names
             #     name_id = findfirst(x -> x == string(t.name), p.parameter_names)
             #     return ASTNode{T}(nothing, PARAMETER, nothing, name_id, nothing, nothing, )
@@ -359,12 +362,22 @@ end
 ########################################################
 # Front facing APIinc
 ########################################################
+"""
+$(TYPEDEF)
+"""
 abstract type AbstractExpressionFunction{T, N, D} <: Function end
 
+"""
+$(TYPEDFIELDS)
+$(TYPEDEF)
+"""
 struct ScalarExpressionFunction{T <: Number} <: AbstractExpressionFunction{T, Node{T, DEFAULT_MAX_DEGREE}, ntuple_type}
     expr::Expression{T, Node{T, DEFAULT_MAX_DEGREE}, ntuple_type}
     num_vars::Int
 
+    """
+    $(TYPEDSIGNATURES)
+    """
     function ScalarExpressionFunction{T}(string::String, var_names::Vector{String}) where T <: Number
         p = Parser{T}(string, var_names)
         # params = _find_parameters(p)
@@ -402,6 +415,11 @@ function (f::ScalarExpressionFunction)(X::SVector{ND, T}, t::T) where {ND, T <:
     return f.expr(vars)[1]
 end
 
+"""
+$(METHODLIST)
+$(TYPEDFIELDS)
+$(TYPEDEF)
+"""
 struct VectorExpressionFunction{N, T <: Number} <: AbstractExpressionFunction{T, Node{T, DEFAULT_MAX_DEGREE}, ntuple_type}
     exprs::SVector{N, ScalarExpressionFunction{T}}
     num_vars::Int
diff --git a/src/Formulations.jl b/src/Formulations.jl
index f692238..4685533 100644
--- a/src/Formulations.jl
+++ b/src/Formulations.jl
@@ -7,6 +7,9 @@ abstract type AbstractElementFormulation{ND, NF} end
 $(TYPEDSIGNATURES)
 """
 num_fields(::AbstractElementFormulation{ND, NF}) where {ND, NF} = NF
+"""
+$(TYPEDSIGNATURES)
+"""
 num_dimensions(::AbstractElementFormulation{ND, NF}) where {ND, NF} = ND
 
 function discrete_gradient end # eventually to be deprecated
diff --git a/src/FunctionSpaces.jl b/src/FunctionSpaces.jl
index 1f9239a..782a9ae 100644
--- a/src/FunctionSpaces.jl
+++ b/src/FunctionSpaces.jl
@@ -104,7 +104,7 @@ default code path that sets up ref fes as a namedtuple
 function _setup_ref_fes(
   mesh::AbstractMesh, 
   interp_type, p_degree,
-  q_type::Type{<:AbstractQuadratureType}, q_degree
+  q_type::Type{<:ReferenceFiniteElements.AbstractQuadratureType}, q_degree
 )
   block_names = mesh.element_block_names
   ref_fes = ReferenceFE[]
@@ -187,7 +187,7 @@ function FunctionSpace(
   is_juliac_safe::Bool = false,
   p_degree::Union{Int, Nothing} = nothing,
   q_degree::Union{Int, Nothing} = nothing,
-) where {IT, QT <: AbstractQuadratureType}
+) where {IT, QT <: ReferenceFiniteElements.AbstractQuadratureType}
   return FunctionSpace{is_juliac_safe}(mesh, field_type, interp_type, QT; p_degree = p_degree, q_degree = q_degree)
 end
 
@@ -196,7 +196,7 @@ function FunctionSpace{is_juliac_safe}(
   q_type::Type{QT} = GaussLobattoLegendre;
   p_degree::Union{Int, Nothing} = nothing,
   q_degree::Union{Int, Nothing} = nothing,
-) where {is_juliac_safe, QT <: AbstractQuadratureType}
+) where {is_juliac_safe, QT <: ReferenceFiniteElements.AbstractQuadratureType}
   # TODO move to some common function so we can use it across
   # all constructors
   if p_degree !== nothing
@@ -238,7 +238,7 @@ function FunctionSpace{is_juliac_safe}(
   q_type::Type{QT} = GaussLobattoLegendre;
   p_degree = nothing,
   q_degree = nothing
-) where {is_juliac_safe, QT <: AbstractQuadratureType}
+) where {is_juliac_safe, QT <: ReferenceFiniteElements.AbstractQuadratureType}
   if is_juliac_safe
     ref_fes = _setup_juliac_safe_ref_fes(interp_type, q_type)
   else