diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 1d306e1..f7dd871 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -19,7 +19,7 @@ jobs:
       matrix:
         version:
           - '1'
-          - '1.8'
+          - '1.11'
         os:
           - ubuntu-latest
           - macOS-latest
diff --git a/.typos.toml b/.typos.toml
index 8a8d8b3..ac18a55 100644
--- a/.typos.toml
+++ b/.typos.toml
@@ -1,2 +1,5 @@
 [files]
 extend-exclude = ["*.toml"]
+
+[default.extend-words]
+OT = "OT"
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 1f1991f..0a0f2ca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,15 +4,21 @@ authors = ["Knut Andreas Meyer and contributors"]
 version = "0.2.4"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 
 [compat]
-julia = "1.8"
+julia = "1.11"
 
 [extras]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2"
+
+[sources]
+MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"}
+#MaterialModelsTesting = {path = "/Users/knutan/.julia/dev/MaterialModelsTesting"}
 
 [targets]
-test = ["Test", "ForwardDiff"]
+test = ["Test", "FiniteDiff", "MaterialModelsTesting"]
diff --git a/docs/make.jl b/docs/make.jl
index 443ef17..ca062e1 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -19,6 +19,7 @@ makedocs(;
         "Conversion" => "conversion.md",
         "Differentiation" => "differentiation.md",
         "Implementation" => "implementing.md",
+        "Devdocs" => "devdocs.md"
     ],
 )
 
diff --git a/docs/src/conversion.md b/docs/src/conversion.md
index 0f841ee..19a2a45 100644
--- a/docs/src/conversion.md
+++ b/docs/src/conversion.md
@@ -2,25 +2,33 @@
 CurrentModule = MaterialModelsBase
 ```
 # Conversion
-`MaterialModelsBase` defines an interface for converting parameters and state variables
+`MaterialModelsBase` defines an interface for converting parameters, state variables, and other objects
 to and from `AbstractVector`s. This is useful when doing parameter identification, or 
 when interfacing with different languages that require information to be passed as arrays.
 
 These function can be divided into those providing information about the material and state variables, and those doing the actual conversion.
 
 ## Information functions
+The following should be defined for new materials,
 ```@docs
 get_tensorbase
+get_vector_length
+get_vector_eltype
+```
+
+Whereas the following already have default implementations that should work provided that the above are implemented.
+```@docs
 get_num_statevars
-get_statevar_eltype
-get_num_params
-get_params_eltype
 get_num_tensorcomponents
 ```
 
 ## Conversion functions
+The following should be defined for new materials, state variables, and alike
 ```@docs
 tovector!
 fromvector
+```
+whereas the following already have default implementations that should work provided that the above functions are implemented,
+```@docs
 tovector
 ```
diff --git a/docs/src/devdocs.md b/docs/src/devdocs.md
new file mode 100644
index 0000000..db1ab35
--- /dev/null
+++ b/docs/src/devdocs.md
@@ -0,0 +1,8 @@
+# Internal documentation
+The following documentation returns to internals of the package and may change
+at any time.
+
+```@docs
+MaterialModelsBase.stress_controlled_indices
+MaterialModelsBase.strain_controlled_indices
+```
diff --git a/src/MaterialModelsBase.jl b/src/MaterialModelsBase.jl
index 3d4a00b..68fd10b 100644
--- a/src/MaterialModelsBase.jl
+++ b/src/MaterialModelsBase.jl
@@ -1,5 +1,6 @@
 module MaterialModelsBase
 using Tensors, StaticArrays
+using ForwardDiff: ForwardDiff
 import Base: @kwdef
 
 # General interface for <:AbstractMaterial
@@ -24,8 +25,9 @@ export update_stress_state!                                 # For nonzero stress
 # For parameter identification and differentiation of materials
 export tovector, tovector!, fromvector                      # Convert to/from `AbstractVector`s
 export get_num_tensorcomponents, get_num_statevars          # Information about the specific material
-export get_num_params, get_params_eltype                    # 
-export MaterialDerivatives, differentiate_material!         # Differentiation routines
+export get_vector_length, get_vector_eltype                 # Get the length and type when converting objects to vectors
+export MaterialDerivatives, StressStateDerivatives          # Derivative collections
+export differentiate_material!                              # Differentiation routines
 export allocate_differentiation_output                      # 
 
 abstract type AbstractMaterial end
@@ -112,14 +114,16 @@ end
 Return the (default) initial `state::AbstractMaterialState` 
 of the material `m`. 
 
-Defaults to the empty `NoMaterialState()`
+Defaults to the empty `NoMaterialState{T}()`, where 
+`T = get_vector_eltype(m)` (which defaults to `Float64`).
 """
-function initial_material_state(::AbstractMaterial) 
-    return NoMaterialState()
+function initial_material_state(m::AbstractMaterial)
+    T = get_vector_eltype(m)
+    return NoMaterialState{T}()
 end
 
 abstract type AbstractMaterialState end
-struct NoMaterialState <: AbstractMaterialState end
+struct NoMaterialState{T} <: AbstractMaterialState end
 
 # Material cache 
 """
@@ -153,8 +157,8 @@ abstract type AbstractExtraOutput end
 struct NoExtraOutput <: AbstractExtraOutput end
 
 include("vector_conversion.jl")
-include("differentiation.jl")
 include("stressiterations.jl")
+include("differentiation.jl")
 include("ErrorExceptions.jl")
 
 end
\ No newline at end of file
diff --git a/src/differentiation.jl b/src/differentiation.jl
index 9e7db4a..fa3acfe 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -1,19 +1,31 @@
 struct MaterialDerivatives{T}
     dσdϵ::Matrix{T}
-    dσdⁿs::Matrix{T}
+    #dσdⁿs::Matrix{T}
     dσdp::Matrix{T}
     dsdϵ::Matrix{T}
-    dsdⁿs::Matrix{T}
+    #dsdⁿs::Matrix{T}
     dsdp::Matrix{T}
 end
 
+function Base.getproperty(d::MaterialDerivatives, key::Symbol)
+    if key === :dσdⁿs || key === :dsdⁿs
+        error("You are probably assuming MaterialModelsBase v0.2 behavior for differentiation")
+    else
+        @inline getfield(d, key)
+    end
+end
+
 """
     MaterialDerivatives(m::AbstractMaterial)
 
 A struct that saves all derivative information using a `Matrix{T}` for each derivative,
-where `T=get_params_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, 
-`get_num_statevars`, and `get_num_params`. The values should be updated in `differentiate_material!`
-by direct access of the fields, where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current 
+where `T=get_vector_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, 
+`get_num_statevars`, and `get_vector_length`. `m` must support `tovector` and `fromvector`, while 
+the output of `initial_material_state` must support `tovector`, and in addition the element type
+of `tovector(initial_material_state(m))` must respect the element type in `tovector(m)` for any `m`.
+
+The values should be updated in `differentiate_material!` by direct access of the fields, 
+where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current 
 and old state variables, and `p` the material parameter vector.
 
 * `dσdϵ`
@@ -25,17 +37,16 @@ and old state variables, and `p` the material parameter vector.
 
 """
 function MaterialDerivatives(m::AbstractMaterial)
-    T = get_params_eltype(m)
+    T = get_vector_eltype(m)
     n_tensor = get_num_tensorcomponents(m)
     n_state = get_num_statevars(m)
-    n_params = get_num_params(m)
+    n_params = get_vector_length(m)
+    dsdp = ForwardDiff.jacobian(p -> tovector(initial_material_state(fromvector(p, m))), tovector(m))
     return MaterialDerivatives(
         zeros(T, n_tensor, n_tensor),  # dσdϵ
-        zeros(T, n_tensor, n_state),   # dσdⁿs
         zeros(T, n_tensor, n_params),  # dσdp
         zeros(T, n_state, n_tensor),   # dsdϵ
-        zeros(T, n_state, n_state),    # dsdⁿs
-        zeros(T, n_state, n_params)    # dsdp
+        dsdp
         )
 end
 
@@ -65,5 +76,108 @@ allocate_differentiation_output(::AbstractMaterial) = NoExtraOutput()
 
 Calculate the derivatives and save them in `diff`, see
 [`MaterialDerivatives`](@ref) for a description of the fields in `diff`.
+
+    differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...)
+
+For material models implementing `material_response(m, args...)` and `differentiate_material!(::MaterialDerivatives, m, args...)`,
+this method will work automatically by
+1) Calling `σ, dσdϵ, state = material_response(stress_state, m, args...)` (except that `dσdϵ::FourthOrderTensor{dim = 3}` is extracted)
+2) Calling `differentiate_material!(ssd.mderiv::MaterialDerivatives, m, args..., dσdϵ::FourthOrderTensor{3})`
+3) Updating `ssd` according to the constraints imposed by the `stress_state`.
+
+For material models that directly implement `material_response(stress_state, m, args...)`, this function should be overloaded directly
+to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the 
+(total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated. 
+The exact interface of the `StressStateDerivatives` datastructure is still not fixed, and may change in the future.
 """
-function differentiate_material! end
\ No newline at end of file
+function differentiate_material! end
+
+struct StressStateDerivatives{T}
+    mderiv::MaterialDerivatives{T}
+    dϵdp::Matrix{T}
+    dσdp::Matrix{T}
+    ϵindex::SMatrix{3, 3, Int} # To allow indexing by (i, j) into only
+    σindex::SMatrix{3, 3, Int} # saved values to avoid storing unused rows.
+    # TODO: Reduce the dimensions, for now all entries (even those that are zero) are stored.
+end
+
+function StressStateDerivatives(::AbstractStressState, m::AbstractMaterial)
+    mderiv = MaterialDerivatives(m)
+    np = get_vector_length(m)
+    nt = get_num_tensorcomponents(m) # Should be changed to only save non-controlled entries
+    dϵdp = zeros(nt, np)
+    dσdp = zeros(nt, np)
+    # Should be changed to only save non-controlled entries
+    vo = Tensors.DEFAULT_VOIGT_ORDER[3]
+    if nt == 6
+        index = SMatrix{3, 3, Int}(min(vo[i, j], vo[j, i]) for i in 1:3, j in 1:3)
+    else
+        index = SMatrix{3, 3, Int}(vo)
+    end
+    return StressStateDerivatives(mderiv, dϵdp, dσdp, index, index)
+end
+
+function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where {N}
+    σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response(stress_state, m, ϵ, args...)
+    differentiate_material!(ssd.mderiv, m, ϵ_full, args..., dσdϵ_full)
+    
+    if isa(stress_state, NoIterationState)
+        copy!(ssd.dσdp, ssd.mderiv.dσdp)
+        fill!(ssd.dϵdp, 0)
+    else
+        sc = stress_controlled_indices(stress_state, ϵ)
+        ec = strain_controlled_indices(stress_state, ϵ)
+        dσᶠdϵᶠ_inv = inv(get_unknowns(stress_state, dσdϵ_full)) # f: unknown strain components solved for during stress iterations
+        ssd.dϵdp[sc, :] .= .-dσᶠdϵᶠ_inv * ssd.mderiv.dσdp[sc, :]
+        ssd.dσdp[ec, :] .= ssd.mderiv.dσdp[ec, :] .+ ssd.mderiv.dσdϵ[ec, sc] * ssd.dϵdp[sc, :]
+        ssd.mderiv.dsdp .+= ssd.mderiv.dsdϵ[:, sc] * ssd.dϵdp[sc, :]
+    end
+    return reduce_tensordim(stress_state, σ_full), reduce_stiffness(stress_state, dσdϵ_full), state, ϵ_full
+end
+
+"""
+    stress_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
+
+Get the `N` indices that are stress-controlled in `stress_state`. The tensor input is used to 
+determine if a symmetric or full tensor is used. 
+"""
+function stress_controlled_indices end
+
+"""
+    strain_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
+
+Get the `N` indices that are strain-controlled in `stress_state`. The tensor input is used to 
+determine if a symmetric or full tensor is used. 
+"""
+function strain_controlled_indices end
+
+# NoIterationState
+stress_controlled_indices(::NoIterationState, ::AbstractTensor) = SVector{0,Int}()
+strain_controlled_indices(::NoIterationState, ::SymmetricTensor) = @SVector([1, 2, 3, 4, 5, 6])
+strain_controlled_indices(::NoIterationState, ::Tensor) = @SVector([1, 2, 3, 4, 5, 6, 7, 8, 9])
+
+# UniaxialStress
+stress_controlled_indices(::UniaxialStress, ::SymmetricTensor) = @SVector([2, 3, 4, 5, 6])
+stress_controlled_indices(::UniaxialStress, ::Tensor) = @SVector([2, 3, 4, 5, 6, 7, 8, 9])
+strain_controlled_indices(::UniaxialStress, ::AbstractTensor) = @SVector([1])
+
+# UniaxialNormalStress
+stress_controlled_indices(::UniaxialNormalStress, ::AbstractTensor) = @SVector([2,3])
+strain_controlled_indices(::UniaxialNormalStress, ::SymmetricTensor) = @SVector([1, 4, 5, 6])
+strain_controlled_indices(::UniaxialNormalStress, ::Tensor) = @SVector([1, 4, 5, 6, 7, 8, 9])
+
+# PlaneStress 12 -> 6, 21 -> 9
+stress_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([3, 4, 5])
+stress_controlled_indices(::PlaneStress, ::Tensor) = @SVector([3, 4, 5, 7, 8])
+strain_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([1, 2, 6])
+strain_controlled_indices(::PlaneStress, ::Tensor) = @SVector([1, 2, 6, 9])
+
+# GeneralStressState
+stress_controlled_indices(ss::GeneralStressState{Nσ}, ::AbstractTensor) where Nσ = controlled_indices_from_tensor(ss.σ_ctrl, true, Val(Nσ))
+function strain_controlled_indices(ss::GeneralStressState{Nσ,TT}, ::AbstractTensor) where {Nσ,TT}
+    N = Tensors.n_components(Tensors.get_base(TT)) - Nσ
+    return controlled_indices_from_tensor(ss.σ_ctrl, false, Val(N))
+end
+function controlled_indices_from_tensor(ctrl::AbstractTensor, return_if, ::Val{N}) where N
+    return SVector{N}(i for (i, v) in pairs(tovoigt(SVector, ctrl)) if v == return_if)
+end
\ No newline at end of file
diff --git a/src/stressiterations.jl b/src/stressiterations.jl
index e7b2c9f..dad4a2c 100644
--- a/src/stressiterations.jl
+++ b/src/stressiterations.jl
@@ -25,7 +25,8 @@ See also [`ReducedStressState`](@ref).
 material_response(::AbstractStressState, ::AbstractMaterial, args...)
 
 @inline function material_response(stress_state::AbstractStressState, m::AbstractMaterial, args::Vararg{Any,N}) where N
-    return reduced_material_response(stress_state, m, args...)
+    σ, dσdϵ, state, ϵ_full = stress_state_material_response(stress_state, m, args...)
+    return reduce_tensordim(stress_state, σ), reduce_stiffness(stress_state, dσdϵ), state, ϵ_full
 end
 
 update_stress_state!(::AbstractStressState, σ) = nothing
@@ -37,8 +38,8 @@ Creates a subtype of `AbstractMaterial` that wraps a stress state and a material
 calls to `material_response(w::ReducedStressState, args...)` gives the same result as 
 `material_response(s, m, args...)`. 
 Calls to `initial_material_state`, `allocate_material_cache`, 
-`get_num_tensorcomponents`, `get_num_statevars`, `get_num_params`, 
-`get_params_eltype`, `tovector!`, `tovector`, 
+`get_num_tensorcomponents`, `get_num_statevars`, `get_vector_length`, 
+`get_vector_eltype`, `tovector!`, `tovector`, 
 and `allocate_differentiation_output` are forwarded with `m` as the argument. 
 `fromvector` returns `ReducedStressState` and is supported as well.
 """
@@ -47,8 +48,8 @@ struct ReducedStressState{S<:AbstractStressState,M<:AbstractMaterial} <: Abstrac
     material::M
 end
 for op in ( :initial_material_state, :allocate_material_cache, 
-            :get_num_tensorcomponents, :get_num_statevars, :get_num_params, 
-            :get_params_eltype, :allocate_differentiation_output)
+            :get_num_tensorcomponents, :get_num_statevars, :get_vector_length, 
+            :get_vector_eltype, :allocate_differentiation_output)
     @eval @inline $op(rss::ReducedStressState) = $op(rss.material)
 end
 function tovector!(v::AbstractVector, rss::ReducedStressState)
@@ -168,7 +169,7 @@ get_tolerance(ss::UniaxialNormalStress) = get_tolerance(ss.settings)
 get_max_iter(ss::UniaxialNormalStress) = get_max_iter(ss.settings)
 
 """
-    GeneralStressState(σ_ctrl::AbstractTensor{2,3,Bool}, σ::AbstractTensor{2,3,Bool}; kwargs...)
+    GeneralStressState(σ_ctrl::AbstractTensor{2,3,Bool}, σ::AbstractTensor{2,3}; kwargs...)
 
 Construct a general stress state controlled by `σ_ctrl` whose component is `true` if that 
 component is stress-controlled and `false` if it is strain-controlled. If stress-controlled,
@@ -252,15 +253,15 @@ reduce_tensordim(::Val{dim}, a::SymmetricTensor{2}) where dim = SymmetricTensor{
 reduce_tensordim(::Val{dim}, A::SymmetricTensor{4}) where dim = SymmetricTensor{4,dim}((i,j,k,l)->A[i,j,k,l])
 
 
-function reduced_material_response(stress_state::NoIterationState, 
+function stress_state_material_response(stress_state::NoIterationState, 
         m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N
 
     ϵ_full = expand_tensordim(stress_state, ϵ)
     σ, dσdϵ, new_state = material_response(m, ϵ_full, args...)
-    return reduce_tensordim(stress_state, σ), reduce_tensordim(stress_state, dσdϵ), new_state, ϵ_full
+    return σ, dσdϵ, new_state, ϵ_full
 end
 
-function reduced_material_response(stress_state::IterationState, 
+function stress_state_material_response(stress_state::IterationState, 
         m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N
 
     # Newton options
@@ -274,8 +275,7 @@ function reduced_material_response(stress_state::IterationState,
         σ_mandel = get_unknowns(stress_state, σ_full)
 
         if norm(σ_mandel) < tol
-            dσdϵ = reduce_stiffness(stress_state, dσdϵ_full)
-            return reduce_tensordim(stress_state, σ_full), dσdϵ, new_state, ϵ_full
+            return σ_full, dσdϵ_full, new_state, ϵ_full
         end
 
         dσdϵ_mandel = get_unknowns(stress_state, dσdϵ_full)
@@ -284,14 +284,22 @@ function reduced_material_response(stress_state::IterationState,
     throw(NoStressConvergence("Stress iterations with the NewtonSolver did not converge"))
 end
 
-reduce_stiffness(::State3D, dσdϵ_full::AbstractTensor{4,3}) = dσdϵ_full
+reduce_stiffness(ss::State3D, dσdϵ_full) = reduce_stiffness(Val(3), ss, dσdϵ_full)
+reduce_stiffness(ss::State2D, dσdϵ_full) = reduce_stiffness(Val(2), ss, dσdϵ_full)
+reduce_stiffness(ss::State1D, dσdϵ_full) = reduce_stiffness(Val(1), ss, dσdϵ_full)
 
-function reduce_stiffness(stress_state, dσdϵ_full::AbstractTensor{4,3})
+reduce_stiffness(::Val{3}, ::State3D, dσdϵ_full::AbstractTensor{4,3}) = dσdϵ_full
+
+function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::IterationState, dσdϵ_full::AbstractTensor{4,3})
     ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ = extract_substiffnesses(stress_state, dσdϵ_full)
     dσᶜdϵᶜ = ∂σᶜ∂ϵᶜ - ∂σᶜ∂ϵᶠ * (∂σᶠ∂ϵᶠ \ ∂σᶠ∂ϵᶜ)
     return convert_stiffness(dσᶜdϵᶜ, stress_state, dσdϵ_full)
 end
 
+function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::NoIterationState, dσdϵ_full::AbstractTensor{4,3})
+    return reduce_tensordim(stress_state, dσdϵ_full)
+end
+
 convert_stiffness(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,1}, dσᶜdϵᶜ)
 convert_stiffness(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::Tensor) = frommandel(Tensor{4,1}, dσᶜdϵᶜ)
 convert_stiffness(dσᶜdϵᶜ::SMatrix{3,3}, ::State2D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,2}, dσᶜdϵᶜ)
@@ -404,24 +412,31 @@ end
 
 # GeneralStressState
 function get_full_tensor(state::GeneralStressState{Nσ}, ::TT, v::SVector{Nσ,T}) where {Nσ,T,TT}
-    shear_factor = T(1/√2)
+    TB = Tensors.get_base(TT)
+    shear_factor = if TB == SymmetricTensor{2,3}
+        T(1/√2)
+    elseif TB == Tensor{2,3}
+        one(T)
+    else
+        error("GeneralStressState expects full dimension of strain input")
+    end
     s(i,j) = i==j ? one(T) : shear_factor
     f(i,j) = state.σ_ctrl[i,j] ? v[state.σ_minds[i,j]]*s(i,j) : zero(T)
-    return Tensors.get_base(TT)(f)
+    return TB(f)
 end
 
 function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{2,3,T}) where {Nσ, T}
-    shear_factor = T(√2)
+    shear_factor = a isa SymmetricTensor ? T(√2) : one(T)
     s(i,j) = i==j ? one(T) : shear_factor
     f(c) = ((i,j) = c; s(i,j)*(a[i,j]-state.σ[i,j]))
-    return SVector{Nσ}((f(c) for c in state.σm_inds))
+    return SVector{Nσ,T}((f(c) for c in state.σm_inds))
 end
 
 function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{4,3,T}) where {Nσ,T}
-    shear_factor = T(√2)
+    shear_factor = a isa SymmetricTensor ? T(√2) : one(T)
     s(i,j) = i==j ? one(T) : shear_factor
     f(c1,c2) = ((i,j) = c1; (k,l) = c2; a[i,j,k,l]*s(i,j)*s(k,l))
-    return SMatrix{Nσ,Nσ}((f(c1,c2) for c1 in state.σm_inds, c2 in state.σm_inds))
+    return SMatrix{Nσ,Nσ,T}((f(c1,c2) for c1 in state.σm_inds, c2 in state.σm_inds))
 end
 
 # Extract each part of the stiffness tensor as SMatrix
diff --git a/src/vector_conversion.jl b/src/vector_conversion.jl
index 1160b50..ff1f359 100644
--- a/src/vector_conversion.jl
+++ b/src/vector_conversion.jl
@@ -1,89 +1,115 @@
-
 """
-    get_tensorbase(m::AbstractMaterial)
+    tovector!(v::AbstractVector, obj; offset = 0)
 
-Return the type of the primary input (strain-like) and associated output (stress-like) 
-to the material model. The default is `SymmetricTensor{2, 3}` for small-strain material 
-models in 3D, but it can be 
-* Small strain material model: Small strain tensor, `ϵ::SymmetricTensor{2, dim}`
-* Large strain material model: Deformation gradient, `F::Tensor{2, dim}`
-* Traction-separation law: Displacement jump, `Δu::Vec{dim}`
+Store the parameters in the object `obj` in `v`. This is typically used with `obj` as 
+an `AbstractMaterial`, an `AbstractMaterialState`, or an `AbstractTensor`. The `offset` input makes it possible 
+to write values starting at an offset location in `v`.
 
-Where `dim = 3` is most common, but any value, 1, 2, or 3, is valid. 
+The implementation for `AbstractTensor`s is included in the `MaterialModelsBase.jl` package, and use the Mandel
+notation for the conversion.
 """
-get_tensorbase(::AbstractMaterial) = SymmetricTensor{2, 3}
+function tovector! end
 
 """
-    get_num_tensorcomponents(::AbstractMaterial)
+    tovector([T::Type{<:AbstractArray}], obj)::T
 
-Returns the number of independent components for the given material. 
-
-!!! note
-    It is not required to implement this function, it is inferred by the 
-    implementation of [`get_tensorbase`](@ref)
+Out-of place version of `tovector!`. Relies on `get_vector_length` and 
+`get_vector_eltype` to be correctly defined, and defaults to `T = Array`.
 
+Experimental support for `T = SArray` is also available, but performance may be suboptimal,
+and can be improved by implementing a custom function for the given type.
 """
-get_num_tensorcomponents(m::AbstractMaterial) = Tensors.n_components(get_tensorbase(m))
+@inline tovector(obj) = tovector(Array, obj)
 
-"""
-    get_num_statevars(s::AbstractMaterialState)
-    get_num_statevars(m::AbstractMaterial)
+function tovector(::Type{<:Array}, obj)
+    T = get_vector_eltype(obj)
+    return tovector!(zeros(T, get_vector_length(obj)), obj)
+end
 
-Return the number of state variables.
-A tensorial state variable should be counted by how many components it has. 
-E.g. if a state consists of one scalar and one symmetric 2nd order tensor in 3d,
-`get_num_statevars` should return 7.
+function tovector(::Type{<:SArray}, obj)
+    T = get_vector_eltype(obj)
+    N = get_vector_length(obj)
+    m = MVector{T, N}()
+    @inline tovector!(m, obj)
+    return SVector{T, N}(m)
+end
 
-It suffices to define for the state, `s`, only, but defining for material, `m`, 
-directly as well can improve performance. 
 """
-function get_num_statevars end
+    fromvector(v::AbstractVector, ::OT; offset = 0)
 
-get_num_statevars(::NoMaterialState) = 0
-get_num_statevars(m::AbstractMaterial) = get_num_statevars(initial_material_state(m))
+Output an object of similar type to `OT`, but with parameters according to `v`. This is typically used with
+an `AbstractMaterial`, an `AbstractMaterialState`, or an `AbstractTensor`. The `offset` input makes it possible 
+to read values starting at an offset location in `v`. 
 
+The implementation for `AbstractTensor`s is included in the `MaterialModelsBase.jl` package, and use the Mandel
+notation for the conversion.
 """
-    get_statevar_eltype(s::AbstractMaterialState)
+function fromvector end
 
-Get the type used to store each scalar component of the material state variables,
-defaults to `Float64`.
 """
-get_statevar_eltype(::AbstractMaterialState) = Float64
+    get_vector_length(obj)
 
+Return the length of the vector representation of `obj` when using 
+`tovector` or `tovector!`. The default implementation of `tovector` relies 
+on this function being defined.
 """
-    get_num_params(m::AbstractMaterial)
+function get_vector_length end
 
-Return the number of material parameters in `m`. No default value implemented. 
 """
-function get_num_params end
+    get_vector_eltype(obj)
 
+Return the element type of the vector representation of `obj` when using 
+`tovector`, i.e. with `T = get_vector_eltype(obj)` and `N = get_vector_length(obj)`,
+we have `typeof(obj) == typeof(fromvector(zeros(T, N), obj))`.
 """
-    get_params_eltype(m::AbstractMaterial)
+function get_vector_eltype end
 
-Return the number type for the scalar material parameters, defaults to `Float64`
 """
-get_params_eltype(::AbstractMaterial) = Float64
+    get_tensorbase(m::AbstractMaterial)
 
-# Conversion functions
+Return the type of the primary input (strain-like) and associated output (stress-like) 
+to the material model. The default is `SymmetricTensor{2, 3}` for small-strain material 
+models in 3D, but it can be 
+* Small strain material model: Small strain tensor, `ϵ::SymmetricTensor{2, dim}`
+* Large strain material model: Deformation gradient, `F::Tensor{2, dim}`
+* Traction-separation law: Displacement jump, `Δu::Vec{dim}`
+
+Where `dim = 3` is most common, but any value, 1, 2, or 3, is valid. 
 """
-    tovector!(v::AbstractVector, m::AbstractMaterial; offset = 0)
+get_tensorbase(::AbstractMaterial) = SymmetricTensor{2, 3}
 
-Put the material parameters of `m` into the vector `v`. 
-This is typically used when the parameters should be fitted.
+"""
+    get_num_tensorcomponents(::AbstractMaterial)
 
-    tovector!(v::AbstractVector, s::AbstractMaterialState; offset = 0)
+Returns the number of independent components for the given material. 
 
-Put the state variables in `s` into the vector `v`.
-This is typically used when differentiating the material 
-wrt. the the old state variables.
+!!! note
+    It is not required to implement this function, it is inferred by the 
+    implementation of [`get_tensorbase`](@ref)
 
-    tovector!(v::AbstractVector, a::Union{SecondOrderTensor, FourthOrderTensor}; offset = 0)
+"""
+get_num_tensorcomponents(m::AbstractMaterial) = Tensors.n_components(get_tensorbase(m))
 
-Puts the Mandel components of `a` into the vector `v`.
 """
-function tovector! end
+    get_num_statevars(m::AbstractMaterial)
+
+Return the number of scalar values required to store the state of `m`, i.e. `length(tovector(initial_material_state(m)))`.
+The default implementation works provided that `s = initial_material_state(m)` and `get_vector_length(s)` is defined.
+"""
+function get_num_statevars(m::AbstractMaterial)
+    return get_vector_length(initial_material_state(m))
+end
+
+# Implementations for types defined in MaterialModelsBase
+tovector!(v::AbstractVector, ::NoMaterialState; kwargs...) = v
+fromvector(::AbstractVector{T}, ::NoMaterialState; kwargs...) where {T} = NoMaterialState{T}()
+get_vector_length(::NoMaterialState) = 0
+get_vector_eltype(::NoMaterialState{T}) where {T} = T
 
 # Tensors.jl implementation
+get_vector_length(::TT) where {TT <: AbstractTensor} = Tensors.n_components(Tensors.get_base(TT))
+get_vector_eltype(a::AbstractTensor) = eltype(a)
+
 function tovector!(v::AbstractVector, a::SecondOrderTensor; offset = 0)
     return tomandel!(v, a; offset)
 end
@@ -94,20 +120,6 @@ function tovector!(v::AbstractVector, a::Union{Tensor{4, <:Any, <:Any, M}, Symme
     return v
 end
 
-"""
-    fromvector(v::AbstractVector, ::MT; offset = 0) where {MT<:AbstractMaterial}
-Create a material of type `MT` with the parameters according to `v`
-
-    fromvector(v::AbstractVector, ::ST; offset = 0) where {ST<:AbstractMaterialState}
-Create a material state of type `ST` with the values according to `v`
-
-    fromvector(v, ::TT; offset = 0) where {TT <: Union{SecondOrderTensor, FourthOrderTensor}}
-
-Create a tensor with shape of `TT` with entries from the Mandel components in `v`.
-"""
-function fromvector end
-
-# Tensors.jl implementation
 function fromvector(v::AbstractVector, ::TT; offset = 0) where {TT <: SecondOrderTensor}
     return frommandel(Tensors.get_base(TT), v; offset)
 end
@@ -118,35 +130,3 @@ function fromvector(v::AbstractVector, ::TT; offset = 0) where {TT <: FourthOrde
     N = round(Int, sqrt(M))
     return frommandel(TB, reshape(view(v, offset .+ (1:M)), (N, N)))
 end
-
-"""
-    tovector(m::AbstractMaterial)
-
-Out-of place version of `tovector!`. Relies on `get_num_params` and 
-`get_params_eltype` to be correctly defined
-"""
-function tovector(m::AbstractMaterial)
-    T = get_params_eltype(m)
-    return tovector!(zeros(T, get_num_params(m)), m)
-end
-
-"""
-    tovector(m::AbstractMaterialState)
-
-Out-of place version of `tovector!`. Relies on `get_num_statevars` and 
-`get_statevar_eltype` to be correctly defined
-"""
-function tovector(s::AbstractMaterialState)
-    T = get_statevar_eltype(s)
-    return tovector!(zeros(T, get_num_statevars(s)), s)
-end
-
-tovector(a::Union{SecondOrderTensor, FourthOrderTensor}) = (m = tomandel(a); reshape(m, length(m)))
-
-# Backwards compatibility
-const get_parameter_type = get_params_eltype
-const material2vector! = tovector!
-const material2vector = tovector
-const vector2material = fromvector
-
-export material2vector!, material2vector, vector2material
diff --git a/test/TestMaterials.jl b/test/TestMaterials.jl
deleted file mode 100644
index bff7526..0000000
--- a/test/TestMaterials.jl
+++ /dev/null
@@ -1,87 +0,0 @@
-module TestMaterials
-using MaterialModelsBase
-using Tensors, StaticArrays, ForwardDiff
-
-# Overloaded functions
-import MaterialModelsBase as MMB
-
-# Exported materials
-#export LinearElastic, ViscoElastic
-
-# LinearElastic material 
-struct LinearElastic{T} <: AbstractMaterial 
-    G::T
-    K::T
-end
-function get_stiffness(m::LinearElastic)
-    I2 = one(SymmetricTensor{2,3})
-    Ivol = I2⊗I2
-    Isymdev = minorsymmetric(otimesu(I2,I2) - Ivol/3)
-    return 2*m.G*Isymdev + m.K*Ivol
-end
-
-function MMB.material_response(
-        m::LinearElastic, ϵ::SymmetricTensor{2},
-        old::NoMaterialState, 
-        Δt, cache, extras) # Explicit values added to test the defaults
-    dσdϵ = get_stiffness(m)
-    σ = dσdϵ⊡ϵ
-    return σ, dσdϵ, old
-end
-
-MMB.get_num_params(::LinearElastic) = 2
-function MMB.tovector!(v::Vector, m::LinearElastic; offset = 0)
-    v[1 + offset] = m.G
-    v[2 + offset] = m.K
-    return v
-end
-MMB.fromvector(v::Vector, ::LinearElastic; offset = 0) = LinearElastic(v[1 + offset], v[2 + offset])
-
-function MMB.differentiate_material!(
-    diff::MaterialDerivatives,
-    m::LinearElastic,
-    ϵ::SymmetricTensor,
-    args...)
-    tomandel!(diff.dσdϵ, get_stiffness(m))
-
-    σ_from_p(p::Vector) = tomandel(get_stiffness(fromvector(p, m))⊡ϵ)
-    ForwardDiff.jacobian!(diff.dσdp, σ_from_p, tovector(m))
-end
-
-# ViscoElastic material: Note - not physically reasonable model because 
-# the volumetric response is also viscous
-struct ViscoElastic{T} <: AbstractMaterial
-    E1::LinearElastic{T}
-    E2::LinearElastic{T}
-    η::T
-end
-struct ViscoElasticState{T} <: AbstractMaterialState
-    ϵv::SymmetricTensor{2,3,T,6}
-end
-MMB.initial_material_state(::ViscoElastic) = ViscoElasticState(zero(SymmetricTensor{2,3}))
-
-# To allow easy automatic differentiation
-function get_stress_ϵv(m::ViscoElastic, ϵ::SymmetricTensor{2,3}, old::ViscoElasticState, Δt)
-    E1 = get_stiffness(m.E1)
-    E2 = get_stiffness(m.E2)
-    ϵv = inv(Δt*E2+one(E2)*m.η)⊡(Δt*E2⊡ϵ + m.η*old.ϵv)
-    σ = E1⊡ϵ + m.η*(ϵv-old.ϵv)/Δt
-    return σ, ϵv
-end
-
-function MMB.material_response(
-        m::ViscoElastic, ϵ::SymmetricTensor{2},
-        old::ViscoElasticState, Δt, args...)
-    if Δt == zero(Δt)
-        if norm(ϵ) ≈ 0
-            return zero(ϵ), NaN*zero(get_stiffness(m.E1)), old
-        else
-            throw(NoLocalConvergence("Δt=0 is not allowed for non-zero strain"))
-        end
-    end
-    σ, ϵv = get_stress_ϵv(m,ϵ,old,Δt)
-    dσdϵ = gradient(ϵ_->get_stress_ϵv(m,ϵ_,old,Δt)[1], ϵ)
-    return σ, dσdϵ, ViscoElasticState(ϵv)
-end
-
-end
\ No newline at end of file
diff --git a/test/differentiation.jl b/test/differentiation.jl
index 70ba2f6..4543c5e 100644
--- a/test/differentiation.jl
+++ b/test/differentiation.jl
@@ -1,18 +1,93 @@
+@testset "differentiation_utils" begin
+    for TT in (Tensor, SymmetricTensor)
+        ϵmock = TT == Tensor ? ones(TT{2,3}) : TT{2,3}((i,j) -> i==j ? 1.0 : 1/√2)
+        for stress_state in (
+                UniaxialStress(), UniaxialNormalStress(), PlaneStress(),
+                GeneralStressState(rand(TT{2, 3, Bool}), zero(TT{2, 3})),
+                )
+            sc = MMB.stress_controlled_indices(stress_state, ϵmock)
+            ec = MMB.strain_controlled_indices(stress_state, ϵmock)
+            @test length(union(Set(sc), Set(ec))) == Tensors.n_components(TT{2, 3})
+            @test tomandel(ϵmock)[sc] ≈ MMB.get_unknowns(stress_state, ϵmock)
+        end
+        for stress_state in (UniaxialStrain(), PlaneStrain(), FullStressState())
+            sc = MMB.stress_controlled_indices(stress_state, ϵmock)
+            ec = MMB.strain_controlled_indices(stress_state, ϵmock)
+            @test length(sc) == 0
+            @test length(ec) == Tensors.n_components(TT{2,3})
+        end
+    end
+end
+
 @testset "differentiation" begin
-    G = rand()
-    K = G + rand()
-    m = LinearElastic(G, K)
-    ϵ = rand(SymmetricTensor{2,3})
-    diff = MaterialDerivatives(m)
-    IxI = one(SymmetricTensor{2,3}) ⊗ one(SymmetricTensor{2,3})
-    dσdϵ = 2G*(one(SymmetricTensor{4,3}) - IxI/3) + K*IxI
-    dσdG = 2*(one(SymmetricTensor{4,3}) - IxI/3) ⊡ ϵ
-    dσdK = IxI ⊡ ϵ
+    @testset "Analytical for LinearElastic" begin
+        G = rand()
+        K = G + rand()
+        m = LinearElastic(G, K)
+        ϵ = rand(SymmetricTensor{2,3})
+        diff = MaterialDerivatives(m)
+        IxI = one(SymmetricTensor{2,3}) ⊗ one(SymmetricTensor{2,3})
+        dσdϵ = 2G*(one(SymmetricTensor{4,3}) - IxI/3) + K*IxI
+        dσdG = 2*(one(SymmetricTensor{4,3}) - IxI/3) ⊡ ϵ
+        dσdK = IxI ⊡ ϵ
+
+        old = NoMaterialState{Float64}(); Δt = nothing
+        cache = NoMaterialCache(); extras = allocate_differentiation_output(m)
+        differentiate_material!(diff, m, ϵ, old, Δt, cache, extras, dσdϵ)
+        @test diff.dσdϵ ≈ tomandel(dσdϵ)
+        @test diff.dσdp[:,1] ≈ tomandel(dσdG)
+        @test diff.dσdp[:,2] ≈ tomandel(dσdK)
+    end
 
-    old = NoMaterialState(); Δt = nothing
-    cache = NoMaterialCache(); extras = allocate_differentiation_output(m)
-    differentiate_material!(diff, m, ϵ, old, Δt, cache, extras, dσdϵ)
-    @test diff.dσdϵ ≈ tomandel(dσdϵ)
-    @test diff.dσdp[:,1] ≈ tomandel(dσdG)
-    @test diff.dσdp[:,2] ≈ tomandel(dσdK)
-end
\ No newline at end of file
+    @testset "Numerical (3D material)" begin
+        elastic = LinearElastic(0.52, 0.77)
+        viscoelastic = ViscoElastic(elastic, LinearElastic(0.33, 0.54), 0.36)
+        for m in (elastic, viscoelastic)    
+            @testset "Initial response" begin
+                ϵ = rand(SymmetricTensor{2,3}) * 1e-3
+                Δt = 1e-2
+                test_derivative(m, ϵ, initial_material_state(m), Δt; 
+                    numdiffsettings = (fdtype = Val{:central},),
+                    comparesettings = ()
+                    )
+                diff = MaterialDerivatives(m)
+                copy!(diff.dsdp, rand(size(diff.dsdp)...))
+                test_derivative(m, ϵ, initial_material_state(m), Δt; 
+                    numdiffsettings = (fdtype = Val{:central},),
+                    comparesettings = (),
+                    diff)
+            end
+            @testset "After shear loading" begin
+                ϵ21 = 0.01; num_steps = 10; t_end = 0.01
+                stressfun(p) = runstrain(fromvector(p, m), ϵ21, (2, 1), t_end, num_steps)[1]
+                dσ21_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, tovector(m), Val{:central}; relstep = 1e-6)
+                σv, state, dσ21_dp, diff = runstrain_diff(m, ϵ21, (2, 1), t_end, num_steps)
+                @test σv ≈ stressfun(tovector(m))
+                @test dσ21_dp ≈ dσ21_dp_num
+            end
+            @testset "FullStressState" begin
+                ϵ21 = 0.01; num_steps = 10; t_end = 0.01
+                ss = FullStressState()
+                # Check that we get the same result for runstresstate and runstrain
+                σ_ss, s_ss = runstresstate(ss, m, ϵ21, (2, 1), t_end, num_steps)
+                σ, s = runstrain(m, ϵ21, (2, 1), t_end, num_steps)
+                @test σ_ss ≈ σ
+                @test tovector(s_ss) ≈ tovector(s)
+            end
+            for (stress_state, ij) in (
+                    (UniaxialStress(), (1,1)), (UniaxialStrain(), (1,1)), 
+                    (UniaxialNormalStress(), (1,1)), (UniaxialNormalStress(), (2,1)),
+                    (PlaneStress(), (2, 2)), (PlaneStrain(), (2, 1)),
+                    (GeneralStressState(SymmetricTensor{2,3,Bool}((true, false, false, false, true, true)), rand(SymmetricTensor{2,3})), (2,2))
+                    )
+                @testset "$(nameof(typeof(stress_state))), (i,j) = ($(ij[1]), $(ij[2]))" begin
+                    ϵij = 0.01; num_steps = 10; t_end = 0.01
+                    stressfun(p) = runstresstate(stress_state, fromvector(p, m), ϵij, ij, t_end, num_steps)[1]
+                    dσij_dp_num = FiniteDiff.finite_difference_jacobian(stressfun, tovector(m), Val{:central}; relstep = 1e-6)
+                    σv, state, dσij_dp, diff = runstresstate_diff(stress_state, m, ϵij, ij, t_end, num_steps)
+                    @test dσij_dp ≈ dσij_dp_num
+                end
+            end
+        end
+    end
+end
diff --git a/test/errors.jl b/test/errors.jl
index 9743442..94eadd5 100644
--- a/test/errors.jl
+++ b/test/errors.jl
@@ -10,5 +10,5 @@
     m = LinearElastic(80e3, 150e3)
     ss = PlaneStress(;max_iter = 1, tolerance=1e-100)
     ϵ = rand(SymmetricTensor{2,2})
-    @test_throws NoStressConvergence material_response(ss, m, ϵ, NoMaterialState())
+    @test_throws NoStressConvergence material_response(ss, m, ϵ, NoMaterialState{Float64}())
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 0b5dea9..46e747c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,46 +1,12 @@
 using MaterialModelsBase
 using Test
 using Tensors, StaticArrays
-import Base: @kwdef
+import MaterialModelsBase as MMB
+using FiniteDiff: FiniteDiff
+using MaterialModelsTesting:
+    LinearElastic, ViscoElastic, NeoHooke, test_derivative, obtain_numerical_material_derivative!,
+    runstrain, runstrain_diff, runstresstate, runstresstate_diff
 
-const MMB = MaterialModelsBase
-
-# <temporary, move to TestMaterials>
-@kwdef struct StVenant{T} <: AbstractMaterial
-    G::T
-    K::T
-end
-function MMB.material_response(m::StVenant, F::Tensor{2,3}, old, args...)
-    function free_energy(defgrad)
-        E = (tdot(defgrad) - one(defgrad)) / 2
-        Edev = dev(E)
-        return m.G * Edev ⊡ Edev + m.K * tr(E)^2 / 2
-    end
-    ∂P∂F, P, _ = hessian(free_energy, F, :all)
-    return P, ∂P∂F, old
-end
-
-@kwdef struct NeoHooke{T} <: AbstractMaterial
-    G::T
-    K::T
-end
-function MMB.material_response(m::NeoHooke, F::Tensor{2,3}, old, args...)
-    function free_energy(defgrad)
-        C = tdot(defgrad)
-        detC = det(C)
-        ΨG = (m.G / 2) * (tr(C) / cbrt(detC) - 3)
-        ΨK = (m.K / 2) * (sqrt(detC) - 1)^2
-        return ΨG + ΨK
-    end
-    ∂P∂F, P, _ = hessian(free_energy, F, :all)
-    return P, ∂P∂F, old
-end
-
-MMB.get_tensorbase(::Union{StVenant, NeoHooke}) = Tensor{2,3}
-# </temporary>
-
-include("TestMaterials.jl")
-using .TestMaterials: LinearElastic, ViscoElastic
 include("utils4testing.jl")
 
 include("vector_conversion.jl")
diff --git a/test/stressiterations.jl b/test/stressiterations.jl
index 5169f2f..bc474bd 100644
--- a/test/stressiterations.jl
+++ b/test/stressiterations.jl
@@ -1,18 +1,27 @@
 all_states = (
     FullStressState(), UniaxialStrain(), PlaneStrain(),
     UniaxialStress(), PlaneStress(), UniaxialNormalStress())
-iter_states = ( UniaxialStress(), PlaneStress(), UniaxialNormalStress())
+gss_ctrl = SymmetricTensor{2,3,Bool}((true, true, falses(4)...))
+iter_states = ( UniaxialStress(), PlaneStress(), UniaxialNormalStress(), 
+                (SymmetricTensor = GeneralStressState(gss_ctrl, 0.0 * gss_ctrl),
+                 Tensor = GeneralStressState(Tensor{2,3}(gss_ctrl), 0.0 * Tensor{2,3}(gss_ctrl)))
+                 )
 iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]),
                 ([3,4,5,7,8], [3,4,5]),
-                ([2,3], [2,3]))
-
+                ([2,3], [2,3]),
+                ([1,6,9], [1,6]))
 
 @testset "conversions" begin
     # Test that
     # 1) All conversions are invertible (convert back and fourth)
     # 2) Are compatible with the expected mandel dynamic tensors 
-    for stress_state in all_states
+    for _stress_state in all_states
         for TT in (Tensor, SymmetricTensor)
+            stress_state = if isa(_stress_state, NamedTuple)
+                _stress_state[nameof(TT)]
+            else
+                _stress_state
+            end
             ared = rand(TT{2,_getdim(stress_state)})
             afull = MMB.expand_tensordim(stress_state, ared)
             @test _getdim(afull) == 3    # Full dimensional tensor
@@ -26,8 +35,13 @@ iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]),
         end
     end
 
-    for (state, inds) in zip(iter_states, iter_mandel)
+    for (_stress_state, inds) in zip(iter_states, iter_mandel)
         for (TT, ii) in zip((Tensor, SymmetricTensor), inds)
+            state = if isa(_stress_state, NamedTuple)
+                _stress_state[nameof(TT)]
+            else
+                _stress_state
+            end
             a = rand(TT{2,3})
             A = rand(TT{4,3})
             ax = MMB.get_unknowns(state, a)
@@ -46,8 +60,13 @@ iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]),
 end
 
 @testset "stiffness_calculations" begin
-    for (state, inds) in zip(iter_states, iter_mandel)
+    for (_stress_state, inds) in zip(iter_states, iter_mandel)
         for (TT, ii) in zip((Tensor, SymmetricTensor), inds)
+            state = if isa(_stress_state, NamedTuple)
+                _stress_state[nameof(TT)]
+            else
+                _stress_state
+            end
             dσdϵ = rand(TT{4,3}) + one(TT{4,3})
             dσᶜdϵᶜ = MMB.reduce_stiffness(state, dσdϵ)
             if _getdim(state) < 3   # Otherwise, conversion is direct    
@@ -150,30 +169,29 @@ end
 
     Δϵ = 1e-6   # Small strain to compare with linear case
     rtol = 1e-5 # Relative tolerance to compare with linear case
-    @testset for m in (StVenant(;G, K), NeoHooke(;G, K))
-        old = initial_material_state(m)
-        # UniaxialStress
-        P, dPdF, state, Ffull = material_response(UniaxialStress(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0)
-        @test isapprox(P[1,1], E*Δϵ; rtol)
-        @test isapprox(dPdF[1,1,1,1], E; rtol)
-        @test Ffull[2,2] ≈ Ffull[3,3]
-        @test Ffull[2,2] ≈ 1-ν*Δϵ
-        
-        # UniaxialStrain
-        P, dPdF, state, Ffull = material_response(UniaxialStrain(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0)
-        @test isapprox(P[1,1], (2G+λ)*Δϵ; rtol)
-        @test isapprox(dPdF[1,1,1,1], (2G+λ); rtol)
-        P_full, dPdF_full, _ = material_response(m, Ffull, old, 0.0)
-        @test isapprox(P_full[2,2], λ*Δϵ; rtol)
-
-        # PlaneStress
-        ϵ = Δϵ * rand(SymmetricTensor{2,2})
-        F = one(Tensor{2,2}) + ϵ
-        P, dPdF, state, Ffull = material_response(PlaneStress(), m, F, old, 0.0)
-        Dvoigt = (E/(1-ν^2))*[1 ν 0; ν 1 0; 0 0 (1-ν)/2]
-        @test isapprox(tovoigt(symmetric(dPdF)), Dvoigt; rtol)
-        @test isapprox(P, fromvoigt(SymmetricTensor{4,2}, Dvoigt)⊡ϵ; rtol)
-    end
+    m = NeoHooke(;G, K)
+    old = initial_material_state(m)
+    # UniaxialStress
+    P, dPdF, state, Ffull = material_response(UniaxialStress(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0)
+    @test isapprox(P[1,1], E*Δϵ; rtol)
+    @test isapprox(dPdF[1,1,1,1], E; rtol)
+    @test Ffull[2,2] ≈ Ffull[3,3]
+    @test Ffull[2,2] ≈ 1-ν*Δϵ
+    
+    # UniaxialStrain
+    P, dPdF, state, Ffull = material_response(UniaxialStrain(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0)
+    @test isapprox(P[1,1], (2G+λ)*Δϵ; rtol)
+    @test isapprox(dPdF[1,1,1,1], (2G+λ); rtol)
+    P_full, dPdF_full, _ = material_response(m, Ffull, old, 0.0)
+    @test isapprox(P_full[2,2], λ*Δϵ; rtol)
+
+    # PlaneStress
+    ϵ = Δϵ * rand(SymmetricTensor{2,2})
+    F = one(Tensor{2,2}) + ϵ
+    P, dPdF, state, Ffull = material_response(PlaneStress(), m, F, old, 0.0)
+    Dvoigt = (E/(1-ν^2))*[1 ν 0; ν 1 0; 0 0 (1-ν)/2]
+    @test isapprox(tovoigt(symmetric(dPdF)), Dvoigt; rtol)
+    @test isapprox(P, fromvoigt(SymmetricTensor{4,2}, Dvoigt)⊡ϵ; rtol)
 end
 
 @testset "visco_elastic" begin