diff --git a/src/Enzyme.jl b/src/Enzyme.jl
index e6cc8a8..85e6653 100644
--- a/src/Enzyme.jl
+++ b/src/Enzyme.jl
@@ -11,7 +11,6 @@ function assemble_scalar_enzyme_safe!(
   U_old = p.field_old
   _update_for_assembly!(p, dof, Uu)
   conns = fspace.elem_conns
-  # foreach_block(fspace, p) do physics, props, ref_fe, b
   for (b, (
     block_physics, ref_fe
   )) in enumerate(zip(
@@ -20,14 +19,13 @@ function assemble_scalar_enzyme_safe!(
     _assemble_scalar_block_enzyme_safe!(
       KA.CPU(),
       block_view(storage, b),
-      conns.data, conns.offsets[b], 
+      conns, 
       func,
       b,
       block_physics, ref_fe,
       X, t, Δt,
       U, U_old,
-      block_view(p.state_old, b), block_view(p.state_new, b),
-      p.properties,
+      p.state_old, p.state_new, p.properties,
       return_type
     )
   end
@@ -35,49 +33,35 @@ end
 
 function _assemble_scalar_block_enzyme_safe!(
   ::KA.CPU,
-  # field::AbstractField,
   field,
-  conns::Conn, coffset::Int,
+  conns_all,
   func::Function,
   b::Int,
   physics::AbstractPhysics, ref_fe::ReferenceFE,
   X::AbstractField, t::T, dt::T,
   U::Solution, U_old::Solution, 
-  state_old::S, state_new::S, props::AbstractArray,
+  state_old::StateVariableField, state_new::StateVariableField, props::PropertyField,
   return_type::R
 ) where {
   T        <: Number,
-  Conn     <: AbstractArray,
   Solution <: AbstractField,
-  S,       #<: L2QuadratureField
   R        <: AssembledReturnType
 }
 
-  for e in axes(state_old, 3)
+  conns = conns_all.data
+  coffset = conns_all.offsets[b]
+  for e in 1:conns_all.nelems[b]
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old)
     props_el = properties(props, e, b)
-    # val_el = _element_scratch(return_type, ref_fe, U)
 
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)
-      state_old_q = _quadrature_level_state(state_old, q, e)
-      state_new_q = _quadrature_level_state(state_new, q, e)
+      state_old_q = state_variables(state_old, q, e, b)
+      state_new_q = state_variables(state_new, q, e, b)
       val_q = func(physics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el)
-      # val_el = _accumulate_q_value(return_type, field, val_q, val_el, q, e)
       field[1, q, e] = val_q
     end
-    # _assemble_element!(field, val_el, conn, e)
-    
-    # writing inline to avoid atomic call
-    # n_dofs = size(field, 1)
-    # for d in axes(field, 1)
-    #   for n in axes(conn, 1)
-    #     global_id = n_dofs * (conn[n] - 1) + d
-    #     local_id = n_dofs * (n - 1) + d
-    #     field.data[global_id] += val_el[local_id]
-    #   end
-    # end
   end
   return nothing
 end
@@ -118,18 +102,15 @@ function assemble_vector_enzyme_safe!(
     values(p.physics), values(fspace.ref_fes)
   ))
     _assemble_vector_block_enzyme_safe!(
-      # KA.get_backend(storage),
       KA.CPU(),
       storage,
-      conns.data, conns.offsets[b], 
+      conns,
       func,
       b,
       block_physics, ref_fe,
       X, t, Δt,
       U, U_old,
-      block_view(p.state_old, b), block_view(p.state_new, b),
-      p.properties
-      # return_type
+      p.state_old, p.state_new, p.properties
     )
   end
   
@@ -146,25 +127,21 @@ TODO add state variables and physics properties
 """
 function _assemble_vector_block_enzyme_safe!(
   ::KA.CPU,
-  # field::AbstractField,
   field,
-  conns::Conn, coffset::Int,
+  conns_all,
   func::Function,
   b::Int,
   physics::AbstractPhysics, ref_fe::ReferenceFE,
   X::AbstractField, t::T, dt::T,
   U::Solution, U_old::Solution, 
-  state_old::S, state_new::S, props::AbstractArray,
-  # return_type::R
+  state_old::StateVariableField, state_new::StateVariableField, props::PropertyField,
 ) where {
   T        <: Number,
-  Conn     <: AbstractArray,
-  Solution <: AbstractField,
-  S,       #<: L2QuadratureField
-  # R        <: AssembledReturnType
+  Solution <: AbstractField
 }
-
-  for e in axes(state_old, 3)
+  conns = conns_all.data
+  coffset = conns_all.offsets[b]
+  for e in 1:conns_all.nelems[b]
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old)
     
@@ -173,8 +150,8 @@ function _assemble_vector_block_enzyme_safe!(
 
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)
-      state_old_q = _quadrature_level_state(state_old, q, e)
-      state_new_q = _quadrature_level_state(state_new, q, e)
+      state_old_q = state_variables(state_old, q, e, b)
+      state_new_q = state_variables(state_new, q, e, b)
       # val_q = func(physics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el)
       # val_el = _accumulate_q_value(return_type, field, val_q, val_el, q, e)
 
diff --git a/src/Fields.jl b/src/Fields.jl
index 3d2841f..b803f81 100644
--- a/src/Fields.jl
+++ b/src/Fields.jl
@@ -548,11 +548,17 @@ function num_fields(field::PropertyField, b::Int)
 end
 
 function properties(field::PropertyField, e::Int, b::Int)
+    @assert 1 <= b && b <= field.nblocks
     offset = field.offsets[b]
     nfields = num_fields(field, b)
+    # `if/elseif` with no `else` leaves `start` undefined on any third value,
+    # which surfaces as an `UndefVarError` rather than saying what went wrong.
+    # Only two layouts exist, so make the second branch total and assert it.
     if field.isblockconstant[b] == PROPS_CONST
         start = offset
-    elseif field.isblockconstant[b] == PROPS_ELEMS
+    else
+        @assert field.isblockconstant[b] == PROPS_ELEMS
+        @assert 1 <= e && e <= field.nelems[b]
         start = offset + nfields * (e - 1)
     end
     return PropertyFieldView(field.data, start, nfields)
@@ -564,30 +570,31 @@ struct PropertyFieldView{T, D <: AbstractVector{T}} <: AbstractVector{T}
     len::Int
 end
 
-Base.size(v::PropertyFieldView) = (v.len,)
-Base.length(v::PropertyFieldView) = v.len
-Base.IndexStyle(::Type{<:PropertyFieldView}) = IndexLinear()
 Base.@propagate_inbounds function Base.getindex(v::PropertyFieldView, i::Int)
     @boundscheck checkbounds(v, i)
     return @inbounds v.data[v.start + i - 1]
 end
+Base.IndexStyle(::Type{<:PropertyFieldView}) = IndexLinear()
+Base.length(v::PropertyFieldView) = v.len
+Base.size(v::PropertyFieldView) = (v.len,)
 
 ######################################################################################################
 # StateVariableField
 ######################################################################################################
 struct StateVariableField{
     T, # Let it be anything to allow for structs
-    D <: AbstractVector{T}
+    D <: AbstractVector{T},
+    I <: AbstractVector{Int}
 } <: AbstractDiscontinuousField{T, D}
     data::D                    # flat storage (CPU or GPU)
     nblocks::Int
-    nfields::Vector{Int}
-    nepes::Vector{Int} # num nodes, q points, etc.
-    nelems::Vector{Int}
-    offsets::Vector{Int}
+    nfields::I
+    nepes::I # num nodes, q points, etc.
+    nelems::I
+    offsets::I
 
-    function StateVariableField{T, D}(data, nblocks, nfields, nepes, nelems, offsets) where {T, D}
-        new{T, D}(data, nblocks, nfields, nepes, nelems, offsets)
+    function StateVariableField{T, D, I}(data, nblocks, nfields, nepes, nelems, offsets) where {T, D, I}
+        new{T, D, I}(data, nblocks, nfields, nepes, nelems, offsets)
     end
 
     function StateVariableField(arrs::Vector{<:AbstractArray{T, 3}}) where T
@@ -601,7 +608,7 @@ struct StateVariableField{
             offset += nfields[b] * nepes[b] * nelems[b]
         end
         data = mapreduce(vec, vcat, arrs)
-        return StateVariableField{T, typeof(data)}(data, length(nepes), nfields, nepes, nelems, offsets)
+        return StateVariableField{T, typeof(data), typeof(nepes)}(data, length(nepes), nfields, nepes, nelems, offsets)
     end
 
     function StateVariableField(::UndefInitializer, ::Type{T}, nfields::Int, qsizes::Vector{Tuple{Int, Int}}) where T
@@ -621,15 +628,16 @@ struct StateVariableField{
     end
 end
 
-function Adapt.adapt_structure(to, field::StateVariableField{T, D}) where {T, D}
+function Adapt.adapt_structure(to, field::StateVariableField{T, D, I}) where {T, D, I}
     data = adapt(to, field.data)
-    return StateVariableField{T, typeof(data)}(
+    nfields = adapt(to, field.nfields)
+    return StateVariableField{T, typeof(data), typeof(nfields)}(
         data,
         field.nblocks,
-        field.nfields,
-        field.nepes,
-        field.nelems,
-        field.offsets
+        nfields,
+        adapt(to, field.nepes),
+        adapt(to, field.nelems),
+        adapt(to, field.offsets)
     )
 end
 
@@ -648,3 +656,32 @@ end
 function num_fields(field::StateVariableField, b::Int)
     return field.nfields[b]
 end
+
+function state_variables(field::StateVariableField, q::Int, e::Int, b::Int)
+    @assert 1 <= q <= field.nepes[b]
+    @assert 1 <= e <= field.nelems[b]
+    offset  = field.offsets[b]
+    nfields = field.nfields[b]
+    nqs     = field.nepes[b]
+    start   = offset + nfields * (q - 1) + nfields * nqs * (e - 1)
+    return StateVariableFieldView(field.data, start, nfields)
+end
+
+struct StateVariableFieldView{T, D <: AbstractVector{T}} <: AbstractVector{T}
+    data::D
+    start::Int
+    len::Int
+end
+
+Base.@propagate_inbounds function Base.getindex(v::StateVariableFieldView, i::Int)
+    @boundscheck checkbounds(v, i)
+    return @inbounds v.data[v.start + i - 1]
+end
+Base.IndexStyle(::Type{<:StateVariableFieldView}) = IndexLinear()
+Base.length(v::StateVariableFieldView) = v.len
+Base.@propagate_inbounds function Base.setindex!(v::StateVariableFieldView{T, D}, val::T, i::Int) where {T, D}
+    @boundscheck checkbounds(v, i)
+    @inbounds v.data[v.start + i - 1] = val
+    return nothing
+end
+Base.size(v::StateVariableFieldView) = (v.len,)
diff --git a/src/Utils.jl b/src/Utils.jl
index 257e90d..6a8e96c 100644
--- a/src/Utils.jl
+++ b/src/Utils.jl
@@ -267,6 +267,32 @@ end
     quote $(exprs...) end
 end
 
+"""
+takes in connectivity and element id
+"""
+function foreach_element(
+    f, conns, block_id,
+    backend = KA.get_backend(conns.data);
+    max_tasks = Threads.nthreads(),
+    min_elems = 1,
+    prefer_threads::Bool = true,
+    # GPU settings
+    block_size = 256
+)
+    nelems = conns.nelems[block_id]
+    if AK.use_gpu_algorithm(backend, prefer_threads)
+        AK._forindices_gpu(f, 1:nelems, backend; block_size)
+    elseif max_tasks == 1
+        _forindices_serial(f, 1:nelems)
+    else
+        _forindices_threads(
+            f, 1:nelems,
+            max_tasks = max_tasks,
+            min_elems = min_elems
+        )
+    end
+end
+
 #########################################
 # hooks for extensions
 #########################################
diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl
index e6580b1..9a34da8 100644
--- a/src/assemblers/Assemblers.jl
+++ b/src/assemblers/Assemblers.jl
@@ -187,36 +187,6 @@ end
   return x_el, u_el, u_el_old, v_el
 end
 
-"""
-$(TYPEDSIGNATURES)
-"""
-# GPU safe: the property count comes from the physics type, so the slice is a
-# statically sized SVector.  `view(field.data, range)` builds a SubArray whose
-# construction lowers to a dynamic call inside a device kernel.
-@inline function properties(
-  field::PropertyField, ::AbstractPhysics{NF, NP, NS}, e::Int, b::Int
-) where {NF, NP, NS}
-  offset = field.offsets[b]
-  ec = ifelse(field.isblockconstant[b] == PROPS_CONST, 1, e)
-  base = offset + NP * (ec - 1)
-  # return SVector{NP, eltype(field)}(ntuple(i -> field.data[base + i - 1], NP))
-  nprops = num_fields(field, b)
-  # return view(field.data, base:base + NP - 1)
-  # return PropertyFieldView(field.data, base, NP)
-  return PropertyFieldView(field.data, base, nprops)
-end
-
-@inline function _element_level_properties(props::AbstractArray, ::Int)
-  return props
-end
-
-"""
-$(TYPEDSIGNATURES)
-"""
-@inline function _element_level_properties(props::SVector{NP, T}, ::Int) where {NP, T}
-  return props
-end
-
 """
 $(TYPEDSIGNATURES)
 """
@@ -258,13 +228,13 @@ end
   return zeros(SVector{NxNDof, eltype(U)})
 end
 
-"""
-$(TYPEDSIGNATURES)
-"""
-function _quadrature_level_state(state::AbstractArray{<:Number, 3}, q::Int, e::Int)
-  state_q = view(state, :, q, e)
-  return state_q
-end
+# """
+# $(TYPEDSIGNATURES)
+# """
+# function _quadrature_level_state(state::AbstractArray{<:Number, 3}, q::Int, e::Int)
+#   state_q = view(state, :, q, e)
+#   return state_q
+# end
 
 function _sparse_matrix_mass(asm::AbstractAssembler, coo_storage)
   type = _sparse_matrix_type(asm)
@@ -417,30 +387,30 @@ end
 
 function _assemble_block!(
   field,
-  conns::Conn, coffset::Int,
+  conns_all,
   func::Function,
   b::Int, # block index
   physics::AbstractPhysics, ref_fe::ReferenceFE,
   X::AbstractField, t::T, dt::T,
   U::Solution, U_old::Solution, 
-  state_old::S, state_new::S, props::AbstractArray,
+  state_old::StateVariableField, state_new::StateVariableField, props::PropertyField,
   return_type::R
 ) where {
   T        <: Number,
-  Conn     <: AbstractArray,
   Solution <: AbstractField,
-  S,       #<: L2QuadratureField
   R        <: AssembledReturnType
 }
-  fec_foraxes(state_old, 3) do e
+  conns = conns_all.data
+  coffset = conns_all.offsets[b]
+  foreach_element(conns_all, b) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old)
     props_el = properties(props, e, b)
     val_el = _element_scratch(return_type, ref_fe, U)
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)
-      state_old_q = _quadrature_level_state(state_old, q, e)
-      state_new_q = _quadrature_level_state(state_new, q, e)
+      state_old_q = state_variables(state_old, q, e, b)
+      state_new_q = state_variables(state_new, q, e, b)
       val_q = func(physics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el)
       val_el = _accumulate_q_value(return_type, field, val_q, val_el, q, e)
     end
@@ -454,24 +424,23 @@ function _assemble_block!(
   b::Int, # block index
   physics::AbstractPhysics,
   t::T, Δt::T,
-  props::P, state_old::S, state_new::S,
-  conns::Conn, coffset::Int, ref_fe::ReferenceFE,
+  props::PropertyField, state_old::StateVariableField, state_new::StateVariableField,
+  conns_all, ref_fe::ReferenceFE,
   X::AbstractField, U::Solution, U_old::Solution
 ) where {
   T        <: Number,
-  S,
-  P,
-  Conn     <: AbstractArray,
   Solution <: AbstractField
 }
-  fec_foraxes(state_old, 3) do e
+  conns = conns_all.data
+  coffset = conns_all.offsets[b]
+  foreach_element(conns_all, b) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old)
     props_el = properties(props, e, b)
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)
-      state_old_q = _quadrature_level_state(state_old, q, e)
-      state_new_q = _quadrature_level_state(state_new, q, e)
+      state_old_q = state_variables(state_old, q, e, b)
+      state_new_q = state_variables(state_new, q, e, b)
       func!(
         field, e, physics, t, Δt, 
         props_el, state_old_q, state_new_q,
diff --git a/src/assemblers/Diagonal.jl b/src/assemblers/Diagonal.jl
index d3ecbb2..c0b58eb 100644
--- a/src/assemblers/Diagonal.jl
+++ b/src/assemblers/Diagonal.jl
@@ -24,16 +24,11 @@ function assemble_diagonal!(
   )
 end
 
-# function assemble_diagonal!(
-#   assembler, func::F, Uu, p
-# ) where F <: Function
-  # storage = assembler.residual_storage
 function assemble_diagonal!(
   storage, pattern, dof, func, Uu, p;
   use_inplace_methods::Bool = false
 )
   fill!(storage, zero(eltype(storage)))
-  # dof = assembler.dof
   fspace = function_space(dof)
   X = coordinates(p)
   t = current_time(p)
@@ -51,21 +46,19 @@ function assemble_diagonal!(
         b,
         physics,
         t, Δt,
-        p.properties,
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        conns.data, conns.offsets[b], ref_fe, X, U, U_old
+        p.properties, p.state_old, p.state_new,
+        conns, ref_fe, X, U, U_old
       )
     else
       _assemble_block!(
         storage,
-        conns.data, conns.offsets[b],
+        conns,
         func,
         b,
         physics, ref_fe,
         X, t, Δt,
         U, U_old,
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        p.properties,
+        p.state_old, p.state_new, p.properties,
         return_type
       )
     end
diff --git a/src/assemblers/LumpedMass.jl b/src/assemblers/LumpedMass.jl
index 49d38a4..3192ced 100644
--- a/src/assemblers/LumpedMass.jl
+++ b/src/assemblers/LumpedMass.jl
@@ -44,14 +44,13 @@ function assemble_lumped_mass!(
   foreach_block(fspace, p) do physics, ref_fe, b
     _assemble_block!(
       storage,
-      conns.data, conns.offsets[b],
+      conns,
       func,
       b,
       physics, ref_fe,
       X, t, Δt,
       U, U_old,
-      block_view(p.state_old, b), block_view(p.state_new, b),
-      p.properties,
+      p.state_old, p.state_new, p.properties,
       return_type
     )
   end
diff --git a/src/assemblers/Matrix.jl b/src/assemblers/Matrix.jl
index b86e86c..61a0a89 100644
--- a/src/assemblers/Matrix.jl
+++ b/src/assemblers/Matrix.jl
@@ -54,21 +54,20 @@ function assemble_matrix!(
         b,
         physics,
         t, dt,
-        p.properties,
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        conns.data, conns.offsets[b], ref_fe, X, U, U_old
+        p.properties, p.state_old, p.state_new,
+        conns,
+        ref_fe, X, U, U_old
       )
     else
       _assemble_block!(
         block_view(storage, pattern, b),
-        conns.data, conns.offsets[b],
+        conns,
         func,
         b,
         physics, ref_fe,
         X, t, dt,
         U, U_old, 
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        p.properties,
+        p.state_old, p.state_new, p.properties,
         return_type
       )
     end
diff --git a/src/assemblers/MatrixAction.jl b/src/assemblers/MatrixAction.jl
index e45d3cd..ef61ba6 100644
--- a/src/assemblers/MatrixAction.jl
+++ b/src/assemblers/MatrixAction.jl
@@ -36,42 +36,41 @@ function assemble_matrix_free_action!(
   foreach_block(fspace, p) do physics, ref_fe, b
     _assemble_block_matrix_free_action!(
       storage,
-      conns.data, conns.offsets[b],
+      conns,
       func_action,
       b,
       physics, ref_fe,
       X, t, Δt,
       U, U_old, V,
-      block_view(p.state_old, b), block_view(p.state_new, b),
-      p.properties
+      p.state_old, p.state_new, p.properties
     )
   end
 end
 
 function _assemble_block_matrix_free_action!(
   field::AbstractField,
-  conns::Conn, coffset,
+  conns_all,
   func_action::Function,
   b::Int,
   physics::AbstractPhysics, ref_fe::ReferenceFE,
   X::AbstractField, t::T, Δt::T,
   U::Solution, U_old::Solution, V::Solution,
-  state_old::S, state_new::S, props::AbstractArray
+  state_old::StateVariableField, state_new::StateVariableField, props::AbstractArray
 ) where {
   T        <: Number,
-  Conn     <: AbstractArray,
-  Solution <: AbstractField,
-  S
+  Solution <: AbstractField
 }
-  fec_foraxes(state_old, 3) do e
+  conns = conns_all.data
+  coffset = conns_all.offsets[b]
+  foreach_element(conns_all, b) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V)
     props_el = properties(props, e, b)
     Kv_el = _element_scratch(AssembledVector(), ref_fe, U)
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)
-      state_old_q = _quadrature_level_state(state_old, q, e)
-      state_new_q = _quadrature_level_state(state_new, q, e)
+      state_old_q = state_variables(state_old, q, e, b)
+      state_new_q = state_variables(state_new, q, e, b)
       Kv_q = func_action(physics, interps, x_el, t, Δt, u_el, u_el_old, v_el, state_old_q, state_new_q, props_el)
       Kv_el = Kv_el + Kv_q
     end
@@ -129,14 +128,13 @@ function assemble_matrix_free_action_full!(
   foreach_block(fspace, p) do physics, ref_fe, b
     _assemble_block_matrix_free_action!(
       storage,
-      conns.data, conns.offsets[b],
+      conns,
       func_action,
       b,
       physics, ref_fe,
       X, t, Δt,
       U, U_old, V,
-      block_view(p.state_old, b), block_view(p.state_new, b),
-      p.properties
+      p.state_old, p.state_new, p.properties
     )
   end
   # The free-DOF entry point above relies on `p.hvp_scratch_field`'s BC
@@ -193,21 +191,20 @@ function assemble_matrix_action!(
         b,
         physics,
         t, Δt,
-        p.properties,
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        conns.data, conns.offsets[b], ref_fe, X, U, U_old, V
+        p.properties, p.state_old, p.state_new,
+        conns,
+        ref_fe, X, U, U_old, V
       )
     else
       _assemble_block_matrix_action!(
         storage,
-        conns.data, conns.offsets[b],
+        conns,
         func,
         b,
         physics, ref_fe,
         X, t, Δt,
         U, U_old, V,
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        p.properties
+        p.state_old, p.state_new, p.properties
       )
     end
   end
@@ -215,29 +212,28 @@ end
 
 function _assemble_block_matrix_action!(
   field::AbstractField, 
-  conns::Conn, coffset,
+  conns_all,
   func::Function,
   b::Int,
   physics::AbstractPhysics, ref_fe::ReferenceFE,
   X::AbstractField, t::T, Δt::T,
   U::Solution, U_old::Solution, V::Solution,
-  state_old::S, state_new::S, props::AbstractArray
+  state_old::StateVariableField, state_new::StateVariableField, props::AbstractArray
 ) where {
   T        <: Number,
-  Conn     <: AbstractArray,
-  Solution <: AbstractField,
-  S        #<: L2QuadratureField
+  Solution <: AbstractField
 }
-  fec_foraxes(state_old, 3) do e
+  conns = conns_all.data
+  coffset = conns_all.offsets[b]
+  foreach_element(conns_all, b) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V)
-    # props_el = _element_level_properties(props, e)
-    props_el = properties(props, physics, e, b)
+    props_el = properties(props, e, b)
     K_el = _element_scratch(AssembledMatrix(), ref_fe, U)
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)
-      state_old_q = _quadrature_level_state(state_old, q, e)
-      state_new_q = _quadrature_level_state(state_new, q, e)
+      state_old_q = state_variables(state_old, q, e, b)
+      state_new_q = state_variables(state_new, q, e, b)
       K_q = func(physics, interps, x_el, t, Δt, u_el, u_el_old, state_old_q, state_new_q, props_el)
       K_el = K_el + K_q
     end
@@ -253,24 +249,23 @@ function _assemble_block_matrix_action!(
   b::Int,
   physics::AbstractPhysics,
   t::T, Δt::T,
-  props::P, state_old::S, state_new::S,
-  conns::Conn, coffset::Int, ref_fe::ReferenceFE,
+  props::PropertyField, state_old::StateVariableField, state_new::StateVariableField,
+  conns_all, ref_fe::ReferenceFE,
   X::AbstractField, U::Solution, U_old::Solution, V::Solution
 ) where {
   T        <: Number,
-  S,
-  P,
-  Conn     <: AbstractArray,
   Solution <: AbstractField
 }
-  fec_foraxes(state_old, 3) do e
+  conns = conns_all.data
+  coffset = conns_all.offsets[b]
+  foreach_element(conns_all, b) do e
     conn = connectivity(ref_fe, conns, e, coffset)
     x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V)
-    props_el = properties(props, physics, e, b)
+    props_el = properties(props, e, b)
     for q in 1:num_cell_quadrature_points(ref_fe)
       interps = _cell_interpolants(ref_fe, q)
-      state_old_q = _quadrature_level_state(state_old, q, e)
-      state_new_q = _quadrature_level_state(state_new, q, e)
+      state_old_q = state_variables(state_old, q, e, b)
+      state_new_q = state_variables(state_new, q, e, b)
       func!(
         field, e, physics, t, Δt, 
         props_el, state_old_q, state_new_q,
diff --git a/src/assemblers/QuadratureQuantity.jl b/src/assemblers/QuadratureQuantity.jl
index d8b0420..daf0398 100644
--- a/src/assemblers/QuadratureQuantity.jl
+++ b/src/assemblers/QuadratureQuantity.jl
@@ -32,14 +32,13 @@ function assemble_quadrature_quantity!(
   foreach_block(fspace, p) do physics, ref_fe, b
     _assemble_block!(
       block_view(storage, b),
-      conns.data, conns.offsets[b], 
+      conns, 
       func,
       b,
       physics, ref_fe,
       X, t, Δt,
       U, U_old,
-      block_view(p.state_old, b), block_view(p.state_new, b),
-      p.properties,
+      p.state_old, p.state_new, p.properties,
       return_type
     )
   end
@@ -69,16 +68,14 @@ function assemble_quadrature_quantity!(
     values(p.physics), values(fspace.ref_fes)
   ))
     _assemble_block!(
-      # backend,
       block_storage,
-      conns.data, conns.offsets[b], 
+      conns,
       func,
       b,
       block_physics, ref_fe,
       X, t, Δt,
       U, U_old,
-      block_view(p.state_old, b), block_view(p.state_new, b),
-      p.properties,
+      p.state_old, p.state_new, p.properties,
       return_type
     )
   end
diff --git a/src/assemblers/Vector.jl b/src/assemblers/Vector.jl
index 443ec30..0cbebe1 100644
--- a/src/assemblers/Vector.jl
+++ b/src/assemblers/Vector.jl
@@ -37,8 +37,6 @@ function assemble_vector!(
   _update_for_assembly!(p, dof, Uu)
   return_type = AssembledVector()
   conns = fspace.elem_conns
-  # foreach_block(conns, p.physics, p.properties, fspace.ref_fes) do physics, props, ref_fe, b
-  # foreach_block(fspace, p) do physics, props, ref_fe, b
   foreach_block(fspace, p) do physics, ref_fe, b
     # if use_sparse_vector
     #   field = block_view(storage, pattern, b)
@@ -54,22 +52,20 @@ function assemble_vector!(
         b,
         physics,
         t, Δt,
-        # props,
-        p.properties,
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        conns.data, conns.offsets[b], ref_fe, X, U, U_old
+        p.properties, p.state_old, p.state_new,
+        conns,
+        ref_fe, X, U, U_old
       )
     else
       _assemble_block!(
         field,
-        conns.data, conns.offsets[b], 
+        conns,
         func,
         b,
         physics, ref_fe,
         X, t, Δt,
         U, U_old,
-        block_view(p.state_old, b), block_view(p.state_new, b),
-        p.properties,
+        p.state_old, p.state_new, p.properties,
         return_type
       )
     end
diff --git a/test/TestFields.jl b/test/TestFields.jl
index 847d453..5bf4042 100644
--- a/test/TestFields.jl
+++ b/test/TestFields.jl
@@ -316,6 +316,60 @@ end
   @test_throws BoundsError block2[5]
 end
 
+@testitem "Fields - test_state_variable_field_view_indexing" begin
+  # `state_variables(field, q, e, b)` must address the same entry that
+  # `block_view(field, b)[:, q, e]` does -- the flat offset arithmetic is the
+  # only thing standing between a quadrature point and its neighbour's state.
+  a1 = rand(2, 3, 5)   # 2 state vars, 3 quadrature points, 5 elements
+  a2 = rand(4, 2, 7)   # a second block with a different shape
+  field = StateVariableField([a1, a2])
+
+  for (b, a) in enumerate((a1, a2))
+    for e in axes(a, 3), q in axes(a, 2)
+      sv = FiniteElementContainers.state_variables(field, q, e, b)
+      @test length(sv) == size(a, 1)
+      @test all(sv[i] ≈ a[i, q, e] for i in axes(a, 1))
+    end
+  end
+
+  # Same eltype trap as PropertyFieldView: an unparameterized `eltype(D)` in
+  # the supertype silently yields `Any`, which un-isbits anything built from it.
+  sv = FiniteElementContainers.state_variables(field, 1, 1, 1)
+  @test eltype(sv) === Float64
+  @test sv isa AbstractVector{Float64}
+  @test eltype(collect(sv)) === Float64
+  @test Base.IndexStyle(typeof(sv)) === IndexLinear()
+end
+
+@testitem "Fields - test_state_variable_field_view_is_bounds_checked" begin
+  # All blocks share one flat vector, so an unchecked read runs into the
+  # neighbouring quadrature point, element, or block rather than failing.
+  field = StateVariableField([rand(2, 3, 5), rand(4, 2, 7)])
+  sv = FiniteElementContainers.state_variables(field, 1, 1, 1)
+  @test length(sv) == 2
+  @test_throws BoundsError sv[3]
+  @test_throws BoundsError sv[0]
+  @test_throws BoundsError sv[3] = 1.0
+end
+
+@testitem "Fields - test_state_variable_field_view_setindex" begin
+  # `setindex!` is what the constitutive update writes new state through, so
+  # it has to land in the flat storage at the same place `getindex` reads.
+  a = rand(2, 3, 4)
+  field = StateVariableField([copy(a)])
+  sv = FiniteElementContainers.state_variables(field, 2, 3, 1)
+  sv[1] = -1.0
+  sv[2] = -2.0
+  @test FiniteElementContainers.block_view(field, 1)[1, 2, 3] ≈ -1.0
+  @test FiniteElementContainers.block_view(field, 1)[2, 2, 3] ≈ -2.0
+  # Nothing else moved.
+  bv = FiniteElementContainers.block_view(field, 1)
+  for e in axes(a, 3), q in axes(a, 2), i in axes(a, 1)
+    (q, e) == (2, 3) && continue
+    @test bv[i, q, e] ≈ a[i, q, e]
+  end
+end
+
 @testitem "Fields - test_state_variable_field" begin
   a1 = rand(2, 3, 40)
   a2 = rand(3, 4, 10)