From 0cec086266dfadced4ea5d26231965e4fb14a422 Mon Sep 17 00:00:00 2001
From: "Craig M. Hamel" <cmhamel32@gmail.com>
Date: Fri, 24 Apr 2026 20:13:11 -0400
Subject: [PATCH 1/2] a lot more stuff for juliac. Got a full app working!

---
 examples/app/Project.toml             |   4 +
 examples/app/input-file.toml          |  18 +--
 examples/app/{build.jl => make.jl}    |   2 +-
 examples/app/src/MyApp.jl             |  82 ++++++++++---
 src/AppTools.jl                       | 142 +++++++++++++++++-----
 src/Fields.jl                         |   2 -
 src/FiniteElementContainers.jl        |   3 +-
 src/FunctionSpaces.jl                 |  33 ++++-
 src/Parameters.jl                     | 167 ++++++++++++++++++++++----
 src/Physics.jl                        | 111 +++++++++--------
 src/Properties.jl                     | 131 ++++++++++++++++++++
 src/Utils.jl                          |  33 ++++-
 src/assemblers/Assemblers.jl          |  10 +-
 src/bcs/DirichletBCs.jl               |   2 +-
 src/bcs/NeumannBCs.jl                 |  17 +--
 src/bcs/RobinBCs.jl                   |  24 ++--
 src/bcs/Sources.jl                    |  14 +--
 test/TestAppTools.jl                  |   2 +-
 test/TestAssemblers.jl                |   2 +-
 test/TestMesh.jl                      |   8 +-
 test/TestProperties.jl                |  25 ++++
 test/input-file.toml                  |   8 +-
 test/mechanics/TestMechanics.jl       |   2 +-
 test/mechanics/TestMechanicsCommon.jl |  16 +--
 test/poisson/TestPoisson.jl           |   4 +-
 test/runtests.jl                      |   2 +-
 26 files changed, 662 insertions(+), 202 deletions(-)
 rename examples/app/{build.jl => make.jl} (95%)
 create mode 100644 src/Properties.jl
 create mode 100644 test/TestProperties.jl

diff --git a/examples/app/Project.toml b/examples/app/Project.toml
index 0413493..efcdd9a 100644
--- a/examples/app/Project.toml
+++ b/examples/app/Project.toml
@@ -3,5 +3,9 @@ name = "MyApp"
 [deps]
 Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
 FiniteElementContainers = "d08262e4-672f-4e7f-a976-f2cea5767631"
+Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ReferenceFiniteElements = "6dc62d09-f8eb-43fd-9672-074e490a997f"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
diff --git a/examples/app/input-file.toml b/examples/app/input-file.toml
index 9075e6a..b489f18 100644
--- a/examples/app/input-file.toml
+++ b/examples/app/input-file.toml
@@ -1,14 +1,14 @@
 [device]
 backend = "cpu"
 
-[functions.zero]
+[functions.zero_ic]
 type = "constant"
 expression = "0.0"
 variables = ["x", "y"]
 
-[functions.one]
+[functions.zero_bc]
 type = "constant"
-expression = "1.0"
+expression = "0.0"
 variables = ["x", "y", "t"]
 
 [mesh]
@@ -16,11 +16,11 @@ file_path = "poisson.g"
 file_type = "exodus"
 
 [[initial_conditions]]
-function = "zero"
-block = "block_1"
-variable = "u"
+function = "zero_ic"
+blocks = ["block_1"]
+variables = ["u"]
 
 [[boundary_conditions.dirichlet]]
-function = "one"
-side_set = "sset_1"
-variable = "u"
+function = "zero_bc"
+side_sets = ["sset_1", "sset_2", "sset_3", "sset_4"]
+variables = ["u"]
diff --git a/examples/app/build.jl b/examples/app/make.jl
similarity index 95%
rename from examples/app/build.jl
rename to examples/app/make.jl
index 2329e6c..4ed4854 100644
--- a/examples/app/build.jl
+++ b/examples/app/make.jl
@@ -10,7 +10,7 @@ img = ImageRecipe(
     file           = "$src_path/src/MyApp.jl",
     trim_mode      = "safe",
     add_ccallables = false,
-    verbose        = true,
+    verbose        = false,
 )
 
 link = LinkRecipe(
diff --git a/examples/app/src/MyApp.jl b/examples/app/src/MyApp.jl
index b2a4091..eb288ee 100644
--- a/examples/app/src/MyApp.jl
+++ b/examples/app/src/MyApp.jl
@@ -2,46 +2,92 @@ import FiniteElementContainers as FEC
 import FiniteElementContainers.AppTools as AT
 using Exodus
 using FiniteElementContainers
+using Krylov
+using LinearAlgebra
+using TimerOutputs
 
 # include files
 include("Physics.jl")
 
-f(X, _) = 2. * π^2 * sin(π * X[1]) * sin(π * X[2])
-bc_zero_func(_, _) = 0.0
+f(X, _) = 2. * π^2 * sin(2π * X[1]) * sin(2π * X[2])
 
 function app_main(ARGS::Vector{String})
     app = AT.App("MyApp")
     AT.add_cli_arg!(app, "--backend"; default = "cpu")
     sim = AT.setup(app, ARGS)
 
+    #####################################
+    # need to define some types
+    #####################################
+    FT = AT.ExpressionFunction{Float64}
+    SPT = FEC.CSCMatrix()
+
     #####################################
     # setup function space
     #####################################
+    # block_ids = FEC._setup_juliac_safe_block_to_ref_fe_id(sim.mesh)
+    # for n in 1:16
+    #     println(Core.stdout, "Block $n element type $(block_ids[n])")
+    # end
     V = FunctionSpace{true}(sim.mesh, H1Field, Lagrange)
     u = ScalarFunction(V, "u")
     dof = DofManager{false}(u)
-    sp_type = FiniteElementContainers.CSCMatrix
-    asm = SparseMatrixAssembler{sp_type, false, false}(dof)
+    asm = SparseMatrixAssembler{SPT, false, false}(dof)
 
-    physics = Poisson(f)
-    props = create_properties(physics)
+    # setup physics and properties
+    physics = Poisson{typeof(f)}[]
+    props = Vector{Float64}[]
+    for _ in 1:length(sim.mesh.element_conns)
+        temp_physics = Poisson(f)
+        push!(physics, temp_physics)
+        push!(props, create_properties(temp_physics))
+    end
 
-    # U = create_unknowns(asm)
-    U = create_field(asm)
+    times = TimeStepper(0.0, 0.0, 1)
+    p = FEC.TypeStableParameters{FT}(sim.mesh, asm, physics, props, sim.ics, sim.dbcs, times)
+    # U = create_field(asm)
+    # Uu = create_unknowns(asm)
 
-    # setting up ics
-    FT = AT.ExpressionFunction{Float64}
-    ics = InitialConditions{FT}(sim.mesh, dof, sim.ics)
-    update_ic_values!(ics, sim.mesh.nodal_coords)
-    update_field_ics!(U, ics)
+    # # setting up ics
+    # # ics = InitialConditions{FT}(sim.mesh, dof, sim.ics)
+    # X = sim.mesh.nodal_coords
+    # # println(Core.stdout, "X size 1 = $(size(X, 1))")
+    # # println(Core.stdout, "X size 2 = $(size(X, 2))")
+    # # for block_id in 1:1
+    # #     ref_fe = FEC.block_reference_element(V, block_id)
+    # #     println(Core.stdout, ref_fe)
+    # # end
+
+    FEC.initialize!(p)
+    # FEC._update_for_assembly!(p, dof, Uu)
 
-    # setting up dbcs
-    dbcs = DirichletBCs{FT}(sim.mesh, dof, sim.dbcs)
-    update_bc_values!(dbcs, sim.mesh.nodal_coords, 1.0)
-    update_field_dirichlet_bcs!(U, dbcs)
+    # # lsolver = DirectLinearSolver(asm)
+    # # lsolver = x -> IterativeLinearSolver(x, :cg)
+    preconditioner = I
+    timer = TimerOutput()
+    workspace = CgWorkspace(stiffness(asm), residual(asm))
+    ΔUu = create_unknowns(asm)
+    lsolver = IterativeLinearSolver(asm, preconditioner, workspace, timer, ΔUu)
+    nlsolver = NewtonSolver(lsolver)
 
-    # # p = create_parameters(mesh, asm, physics, props)
+    # println(Core.stdout, "Number of blocks = $(length(sim.mesh.element_types))")
+    # # N = length(sim.mesh.element_types)
+    # # # N = Val(N)
+    # # N = Val{N}()
+    # # ref_fe = _get_ref_fe(V, N)
+    # # indices = ntuple(x -> x, N)
+    # # indices = ntuple(x -> x, Val(N))
+    # # s = MyStruct{N}(N)::MyStruct{N}
+    integrator = QuasiStaticIntegrator(nlsolver)
     println(sim.log_file.io, "Setup complete")
+    println(sim.log_file.io, "Solving...")
+    evolve!(integrator, p)
+    println(sim.log_file.io, "Solve complete")
+    println(Core.stdout, maximum(p.field.data))
+    println(Core.stdout, minimum(p.field.data))
+
+    # pp = PostProcessor(sim.mesh, "juliac_output.e", u)
+
 end
 
 function @main(ARGS::Vector{String})
diff --git a/src/AppTools.jl b/src/AppTools.jl
index 82dc361..abffd09 100644
--- a/src/AppTools.jl
+++ b/src/AppTools.jl
@@ -14,6 +14,9 @@ import ..FileMesh
 import ..FunctionSpace
 import ..H1Field
 import ..InitialCondition
+import ..NeumannBC
+import ..RobinBC
+import ..Source
 import ..UnstructuredMesh
 import ..element_blocks
 import ..element_ids
@@ -284,6 +287,9 @@ print_dict(l::LogFile, d::Dict{String, String}) = print_dict(l, d)
 #######################################################
 struct BCSettings{T <: Number}
     dirichlet::Vector{DirichletBC{ExpressionFunction{T}}}
+    neumann::Vector{NeumannBC{ExpressionFunction{T}}}
+    robin::Vector{RobinBC{ExpressionFunction{T}}}
+    source::Vector{Source{ExpressionFunction{T}}}
 end
 
 struct FunctionSettings{T <: Number}
@@ -318,35 +324,102 @@ end
 
 function _parse_boundary_condition_settings(log_file, data, functions::FunctionSettings{T}) where T <: Number
     print_banner(log_file, "Boundary conditions")
-    # if length(data["boundary_conditions"])
     if haskey(data, "boundary_conditions")
         bc_settings = data["boundary_conditions"]::Dict{String, Any}
     else
         bc_settings = Dict{String, Any}()
     end
-    dbc_settings = bc_settings["dirichlet"]::Vector{Any}
     dbcs = DirichletBC{ExpressionFunction{T}}[]
-    for bc in dbc_settings
-        # has_key_check(log_file, bc, "function", "key \"function\" not found for DirichletBC")
-        # has_key_check(log_file, bc, "variable", "key \"variable\" not found for DirichletBC")
-        temp = bc::Dict{String, Any}
-        func = temp["function"]::String
-        func = functions.expr_funcs[func]
-        var = temp["variable"]::String
-        if haskey(temp, "block")
-            block = temp["block"]::String
-            push!(dbcs, DirichletBC(var, func; block_name = block))
-        elseif haskey(temp, "node_set")
-            node_set = temp["node_set"]::String
-            push!(dbcs, DirichletBC(var, func; nodeset_name = node_set))
-        elseif haskey(temp, "side_set")
-            side_set = temp["side_set"]::String
-            push!(dbcs, DirichletBC(var, func; sideset_name = side_set))
-        else
-            error_message(log_file, "Requires block, node_set, or side_set")
+    nbcs = NeumannBC{ExpressionFunction{T}}[]
+    rbcs = RobinBC{ExpressionFunction{T}}[]
+    srcs = Source{ExpressionFunction{T}}[]
+    if haskey(bc_settings, "dirichlet")
+        dbc_settings = bc_settings["dirichlet"]::Vector{Any}
+        for bc in dbc_settings
+            # has_key_check(log_file, bc, "function", "key \"function\" not found for DirichletBC")
+            # has_key_check(log_file, bc, "variable", "key \"variable\" not found for DirichletBC")
+            temp = bc::Dict{String, Any}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            vars = temp["variables"]::Vector{String}
+            if haskey(temp, "blocks")
+                blocks = temp["blocks"]::Vector{String}
+                for block in blocks
+                    for var in vars
+                        push!(dbcs, DirichletBC(var, func; block_name = block))
+                    end
+                end
+            elseif haskey(temp, "node_sets")
+                node_sets = temp["node_sets"]::Vector{String}
+                for node_set in node_sets
+                    for var in vars
+                        push!(dbcs, DirichletBC(var, func; nodeset_name = node_set))
+                    end
+                end
+            elseif haskey(temp, "side_sets")
+                side_sets = temp["side_sets"]::Vector{String}
+                for side_set in side_sets
+                    for var in vars
+                        push!(dbcs, DirichletBC(var, func; sideset_name = side_set))
+                    end
+                end
+            else
+                error_message(log_file, "Requires block, node_set, or side_set")
+            end
+        end
+    end
+
+    if haskey(bc_settings, "neumann")
+        nbc_settings = bc_settings["neumann"]::Vector{Any}
+        for bc in nbc_settings
+            temp = bc::Dict{String, Any}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            sidesets = temp["side_sets"]::Vector{String}
+            vars = temp["variables"]::Vector{String}
+            for side_set in sidesets
+                for var in vars
+                    push!(nbcs, NeumannBC(var, func, side_set))
+                end
+            end
         end
     end
-    return BCSettings(dbcs)
+
+    if haskey(bc_settings, "robin")
+        rbc_settings = bc_settings["robin"]::Vector{Any}
+        for bc in rbc_settings
+            temp = bc::Dict{String, Any}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            sidesets = temp["side_sets"]::Vector{String}
+            vars = temp["variables"]::Vector{String}
+            for side_set in sidesets
+                for var in vars
+                    push!(nbcs, RobinBC(var, func, side_set))
+                end
+            end
+        end
+    end
+
+    if haskey(bc_settings, "source")
+        src_settings = bc_settings["source"]::Vector{Any}
+        for src in src_settings
+            temp = src::Dict{String, Any}
+            blocks = temp["blocks"]::Vector{String}
+            func = temp["function"]::String
+            func = functions.expr_funcs[func]
+            sidesets = temp["side_sets"]::Vector{String}
+            vars = temp["variables"]::Vector{String}
+            for block in blocks
+                for var in vars
+                    push!(nbcs, Source(var, func, block))
+                end
+            end
+        end
+    end
+
+    # nbc_settings = 
+    return BCSettings(dbcs, nbcs, rbcs, srcs)
 end
 
 function _parse_function_space_settings(log_file, data)
@@ -395,11 +468,15 @@ function _parse_initial_condition_settings(log_file, data, functions)
     ics = InitialCondition{ExpressionFunction{Float64}}[]
     for ic in ic_settings
         temp = ic::Dict{String, Any}
-        block = temp["block"]::String
+        blocks = temp["blocks"]::Vector{String}
         func = temp["function"]::String
-        var = temp["variable"]::String
         func = functions.expr_funcs[func]
-        push!(ics, InitialCondition(var, func, block))
+        vars = temp["variables"]::Vector{String}
+        for block in blocks
+            for var in vars
+                push!(ics, InitialCondition(var, func, block))
+            end
+        end
     end
     return ics
 end
@@ -585,6 +662,9 @@ struct Simulation{T <: Number, IO, Mesh}
     ics::Vector{InitialCondition{ExpressionFunction{T}}}
     log_file::LogFile{IO}
     mesh::Mesh
+    nbcs::Vector{NeumannBC{ExpressionFunction{T}}}
+    rbcs::Vector{RobinBC{ExpressionFunction{T}}}
+    srcs::Vector{Source{ExpressionFunction{T}}}
 
     function Simulation(settings::InputSettings, log_file::LogFile{IO}) where IO
         print_banner(log_file, "Mesh")
@@ -603,7 +683,10 @@ struct Simulation{T <: Number, IO, Mesh}
         else
             T = Float64
         end
-        new{T, IO, typeof(mesh)}(settings.bcs.dirichlet, settings.ics, log_file, mesh)
+        new{T, IO, typeof(mesh)}(
+            settings.bcs.dirichlet, settings.ics, log_file, mesh,
+            settings.bcs.neumann, settings.bcs.robin, settings.bcs.source
+        )
     end
 end
 
@@ -699,7 +782,6 @@ function generate_app(
 
         function app_main(ARGS::Vector{String})
             app = AT.App(\"$(name)\")
-            AT.add_cli_arg(app, \"--backend\"; default = \"cpu\")
             sim = AT.setup(app, ARGS)
             println(sim.log_file.io, "Setup complete")
         end
@@ -725,14 +807,14 @@ function generate_app(
         src_path = joinpath(@__DIR__)
         @show build_path
         @show src_path
-        rm(build_path; force=true, recursive=true)
+        rm(build_path; force = true, recursive = true)
 
         img = ImageRecipe(
             output_type    = "--output-exe",
             file           = "\$src_path/src/$name.jl",
             trim_mode      = "$trim_str",
             add_ccallables = false,
-            verbose        = true,
+            verbose        = false,
         )
 
         link = LinkRecipe(
@@ -742,9 +824,9 @@ function generate_app(
 
         bun = BundleRecipe(
             link_recipe = link,
-            output_dir  = build_path # or `nothing` to skip bundling
+            #output_dir  = build_path # or `nothing` to skip bundling
+            output_dir  = nothing
         )
-
         compile_products(img)
         link_products(link)
         bundle_products(bun)
diff --git a/src/Fields.jl b/src/Fields.jl
index 24a82e9..d8a2160 100644
--- a/src/Fields.jl
+++ b/src/Fields.jl
@@ -101,7 +101,6 @@ struct Connectivity{
     nepes::Vector{T}
     nelems::Vector{T}
     offsets::Vector{T}
-# end
 
     function Connectivity(mats::Vector{<:AbstractArray{<:Integer, 2}})
         nblocks = length(mats)
@@ -114,7 +113,6 @@ struct Connectivity{
             offset += nepe * nelem
         end
         data = mapreduce(vec, vcat, mats)
-        # return Connectivity(data, nblocks, nepes, nelems, offsets)
         new{eltype(data), typeof(data)}(data, nblocks, nepes, nelems, offsets)
     end
 
diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl
index 7c37daa..50e93b2 100644
--- a/src/FiniteElementContainers.jl
+++ b/src/FiniteElementContainers.jl
@@ -49,6 +49,7 @@ export H1Field
 export HcurlField
 export HdivField
 export L2Field
+export Properties
 export num_entities
 export num_fields
 
@@ -241,7 +242,6 @@ include("meshes/Meshes.jl")
 include("FunctionSpaces.jl")
 include("Functions.jl")
 include("DofManagers.jl")
-
 # where is best to put this?
 # include("Utils.jl")
 
@@ -251,6 +251,7 @@ include("InitialConditions.jl")
 
 include("Formulations.jl")
 include("Physics.jl")
+include("Properties.jl")
 include("assemblers/Assemblers.jl")
 #
 include("TimeSteppers.jl")
diff --git a/src/FunctionSpaces.jl b/src/FunctionSpaces.jl
index 649f759..b5ac77c 100644
--- a/src/FunctionSpaces.jl
+++ b/src/FunctionSpaces.jl
@@ -60,6 +60,8 @@ const _el_name_to_juliac_safe_id = Dict{String, Int}(
   "TETRA10" => 7
 )
 
+const MAX_BLOCKS = 16
+
 function _setup_block_to_ref_fe_id(mesh::AbstractMesh, is_juliac_safe::Bool)
   if is_juliac_safe
     block_ids = Vector{Int}(undef, 0)
@@ -73,6 +75,17 @@ function _setup_block_to_ref_fe_id(mesh::AbstractMesh, is_juliac_safe::Bool)
   end
 end
 
+function _setup_juliac_safe_block_to_ref_fe_id(mesh::AbstractMesh)
+  names = mesh.element_block_names
+  el_types = map(x -> _el_name_to_juliac_safe_id[mesh.element_types[x]], names)
+  N = length(names)
+  return ntuple(i -> i <= N ? el_types[i] : -1, Val(MAX_BLOCKS))  # replace `i` with your actual block value
+end
+
+function _setup_block_to_ref_fe_id(mesh::AbstractMesh)
+  return 1:length(mesh.element_types) |> collect
+end
+
 # Lagrange elements
 # defaulting to fully integrated elements for now
 function _setup_juliac_safe_ref_fes(::Type{Lagrange}, q_type::Type{QT}) where {QT <: AbstractQuadratureType}
@@ -147,13 +160,15 @@ $(TYPEDFIELDS)
 """
 struct FunctionSpace{
   IsJuliaCSafe,
-  IT <: Integer,
-  IV <: AbstractVector{IT},
+  IT            <: Integer,
+  IV            <: AbstractVector{IT},
+  BTRE,
   Coords,
   RefFEs
 } <: AbstractFunctionSpace
   block_names::Vector{String}
-  block_to_ref_fe_id::Vector{IT} # only used in juliac builds
+  # block_to_ref_fe_id::Vector{IT} # only used in juliac builds
+  block_to_ref_fe_id::BTRE
   coords::Coords
   elem_conns::Connectivity{IT, IV}
   elem_id_maps::Vector{Vector{IT}} # TODO create new type for ID map similar to connectivity
@@ -163,7 +178,8 @@ struct FunctionSpace{
     block_names, block_to_ref_fe_id, coords, conns, elem_id_maps, ref_fes
   ) where is_juliac_safe
     new{
-      is_juliac_safe, eltype(conns.data), typeof(conns.data), typeof(coords), typeof(ref_fes)
+      is_juliac_safe, eltype(conns.data), typeof(conns.data), 
+      typeof(block_to_ref_fe_id), typeof(coords), typeof(ref_fes)
     }(block_names, block_to_ref_fe_id, coords, conns, elem_id_maps, ref_fes)
   end
 end
@@ -207,7 +223,12 @@ function FunctionSpace{is_juliac_safe}(
     conns = Connectivity([val for val in values(mesh.element_conns)])
   end
   elem_id_maps = [val for val in values(mesh.element_id_maps)]
-  block_to_ref_fe_id = _setup_block_to_ref_fe_id(mesh, is_juliac_safe)
+  if is_juliac_safe
+    block_to_ref_fe_id = _setup_juliac_safe_block_to_ref_fe_id(mesh)
+  else
+    # block_to_ref_fe_id = _setup_block_to_ref_fe_id(mesh, is_juliac_afe)
+    block_to_ref_fe_id = _setup_block_to_ref_fe_id(mesh)
+  end
 
   return FunctionSpace{is_juliac_safe}(mesh.element_block_names, block_to_ref_fe_id, coords, conns, elem_id_maps, ref_fes)
 end
@@ -273,7 +294,7 @@ function block_reference_element(fspace::FunctionSpace{false, I, V, C, R}, block
   return fspace.ref_fes[block_id]
 end
 
-function block_reference_element_old(fspace::FunctionSpace{true, I, V, C, R}, block_id::Int) where {I, V, C, R}
+function block_reference_element(fspace::FunctionSpace{true, I, V, C, R}, block_id::Int) where {I, V, C, R}
   ref_fe_id = fspace.block_to_ref_fe_id[block_id]
   return fspace.ref_fes[ref_fe_id]
 end
diff --git a/src/Parameters.jl b/src/Parameters.jl
index 58cf3c8..590a4c1 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -1,3 +1,27 @@
+function _setup_state_variables(fspace, physics)
+  state_old = Array{Float64, 3}[]
+  state_new = Array{Float64, 3}[]
+  for (b, val) in enumerate(values(physics))
+    # create state variables for this block physics
+    NS = num_states(val)
+    NQ, NE = block_quadrature_size(fspace, b)
+
+    state_old_temp = zeros(NS, NQ, NE)
+    state_new_temp = zeros(NS, NQ, NE)
+    for e in 1:NE
+      for q in 1:NQ
+        state_old_temp[:, q, e] = create_initial_state(val)
+        state_new_temp[:, q, e] = create_initial_state(val)
+      end
+    end
+    push!(state_old, state_old_temp)
+    push!(state_new, state_new_temp)
+  end
+  state_old = L2Field(state_old)
+  state_new = L2Field(state_new)
+  return state_old, state_new
+end
+
 """
 $(TYPEDEF)
 $(TYPEDSIGNATURES)
@@ -87,26 +111,27 @@ function Parameters(
   end
 
   # setup state variables
-  state_old = Array{Float64, 3}[]
-  state_new = Array{Float64, 3}[]
-  for (b, val) in enumerate(values(physics))
-    # create state variables for this block physics
-    NS = num_states(val)
-    NQ, NE = block_quadrature_size(fspace, b)
-
-    state_old_temp = zeros(NS, NQ, NE)
-    state_new_temp = zeros(NS, NQ, NE)
-    for e in 1:NE
-      for q in 1:NQ
-        state_old_temp[:, q, e] = create_initial_state(val)
-        state_new_temp[:, q, e] = create_initial_state(val)
-      end
-    end
-    push!(state_old, state_old_temp)
-    push!(state_new, state_new_temp)
-  end
-  state_old = L2Field(state_old)
-  state_new = L2Field(state_new)
+  # state_old = Array{Float64, 3}[]
+  # state_new = Array{Float64, 3}[]
+  # for (b, val) in enumerate(values(physics))
+  #   # create state variables for this block physics
+  #   NS = num_states(val)
+  #   NQ, NE = block_quadrature_size(fspace, b)
+
+  #   state_old_temp = zeros(NS, NQ, NE)
+  #   state_new_temp = zeros(NS, NQ, NE)
+  #   for e in 1:NE
+  #     for q in 1:NQ
+  #       state_old_temp[:, q, e] = create_initial_state(val)
+  #       state_new_temp[:, q, e] = create_initial_state(val)
+  #     end
+  #   end
+  #   push!(state_old, state_old_temp)
+  #   push!(state_new, state_new_temp)
+  # end
+  # state_old = L2Field(state_old)
+  # state_new = L2Field(state_new)
+  state_old, state_new = _setup_state_variables(fspace, physics)
 
   # scratch
   field = create_field(dof)
@@ -182,14 +207,87 @@ function KA.get_backend(p::Parameters)
   return KA.get_backend(p.field)
 end
 
+struct TypeStableParameters{
+  # Funcs  <: AbstractVector,
+  FuncT,
+  IT     <: Integer,
+  RT     <: Number,
+  IV     <: AbstractVector{IT},
+  RV     <: AbstractVector{RT},
+  # RM1    <: AbstractMatrix,
+  # RM2    <: AbstractMatrix,
+  # RM3    <: AbstractMatrix,
+  # RM4    <: AbstractMatrix,
+  Phys,
+  Props,
+  Coords <: AbstractField,
+  Field  <: AbstractField 
+} <: AbstractParameters
+  ics::InitialConditions{Vector{InitialConditionFunction{FuncT}}, IV, RV}
+  dirichlet_bcs::DirichletBCs{Vector{DirichletBCFunction{FuncT, FuncT, FuncT}}, IV, RV}
+  # neumann_bcs::NeumannBCs{NBCFuncs, IT, IV, RM1}
+  # robin_bcs::RobinBCs{RBCFuncs, IT, IV, RM2, RM3}
+  # sources::Sources{SRCFuncs, RM4}
+  times::TimeStepper{RT}
+  physics::Phys
+  properties::Props
+  state_old::L2Field{RT, RV}
+  state_new::L2Field{RT, RV}
+  coords::Coords
+  field::Field
+  field_old::Field
+  # scratch fields
+  hvp_scratch_field::Field
+
+  function TypeStableParameters{F}(mesh, assembler, physics, props, ics, dbcs, times) where F
+    dof = assembler.dof
+    fspace = function_space(dof)
+    ics = InitialConditions{F}(mesh, dof, ics)
+    dbcs = DirichletBCs{F}(mesh, dof, dbcs)
+
+    state_old, state_new = _setup_state_variables(fspace, physics)
+
+    coords = mesh.nodal_coords
+    field = create_field(assembler)
+    field_old = create_field(assembler)
+    hvp_scratch_field = create_field(assembler)
+
+    # update assembler, where should this really live?
+    update_dofs!(assembler, dbcs)
+
+    new{
+      F, Int, Float64, Vector{Int}, Vector{Float64},
+      typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field)
+    }(
+      ics, dbcs, times, 
+      physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field
+    )
+
+    # TODO need to settle on expected format for nbcs, rbcs, and sources...
+    # should we always expect a staticarrays even for a 1 dof problem?
+    # RM1 = Matrix{Float64}
+    # RM2 = Matrix{Float64}
+    # RM3 = Matrix{Float64}
+    # RM4 = Matrix{Float64}
+
+    # new{
+    #   F, Int, Float64, Vector{Int}, Vector{Float64}, RM1, RM2, RM3, RM4,
+    #   typeof(physics), typeof(props), typeof(mesh.nodal_coords), typeof(field)
+    # }(
+    #   ics, dbcs, times, 
+    #   physics, props, state_old, state_new, coords, field, field_old, hvp_scratch_field
+    # )
+  end
+end
+
 function create_parameters(
   mesh, assembler, physics, props;
-  ics=InitialCondition[],
-  dirichlet_bcs=DirichletBC[],
-  neumann_bcs=NeumannBC[],
-  robin_bcs=RobinBC[],
-  sources=Source[],
-  times=nothing
+  ics                  = InitialCondition[],
+  dirichlet_bcs        = DirichletBC[],
+  neumann_bcs          = NeumannBC[],
+  robin_bcs            = RobinBC[],
+  sources              = Source[],
+  times                = nothing
 )
   return Parameters(mesh, assembler, physics, props, ics, dirichlet_bcs, neumann_bcs, robin_bcs, sources, times)
 end
@@ -258,6 +356,23 @@ function update_bc_values!(p::AbstractParameters, assembler)
   return nothing
 end
 
+function update_bc_values!(p::TypeStableParameters, assembler)
+  X = coordinates(p)
+  t = current_time(p)
+  update_bc_values!(p.dirichlet_bcs, X, t)
+  # update_bc_values!(p.neumann_bcs, assembler, X, t)
+  # update_source_values!(p.sources, assembler, X, t)
+
+  # TODO how to handle Robin BCs?
+  # currently assembly methods handle updating the field
+  # in parameters with the current unknown dofs
+  # we need field here to reflect that for the robin bcs
+  # to be correct...
+
+  # order of operations goes
+  return nothing
+end
+
 """
 $(TYPEDSIGNATURES)
 """
diff --git a/src/Physics.jl b/src/Physics.jl
index a4c2e21..e3ebcb4 100644
--- a/src/Physics.jl
+++ b/src/Physics.jl
@@ -1,9 +1,38 @@
 abstract type AbstractPhysics{NF, NP, NS} end
-
 num_fields(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS} = NF
 num_properties(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS} = NP
 num_states(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS} = NS
 
+# physics like methods
+function damping end
+function energy end
+function mass end
+function mass! end
+function mass_action end
+function mass_action! end
+function residual end
+function residual! end
+function stiffness end
+function stiffness! end
+function stiffness_action end
+function stiffness_action! end
+
+# optimization like methods
+function gradient end
+function hessian end
+function value end
+
+# default
+function create_function(fspace, physics::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS}
+  names = field_names(physics)
+  return GeneralFunction(fspace, names)
+end
+
+# default
+function create_initial_state(::AbstractPhysics{NF, NP, 0}) where {NF, NP}
+  return SVector{0, Float64}()
+end
+
 """
 $(TYPEDSIGNATURES)
 """
@@ -12,12 +41,19 @@ function create_properties(physics::AbstractPhysics{NF, NP, NS}) where {NF, NP,
   type $(typeof(physics))!"
 end
 
+# default
 function create_properties(::AbstractPhysics{NF, 0, NS}) where {NF, NS}
   return SVector{0, Float64}()
 end
 
-function create_initial_state(::AbstractPhysics{NF, NP, 0}) where {NF, NP}
-  return SVector{0, Float64}()
+# default
+function field_names(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS}
+  if NF == 0
+    field_names = String[]
+  else
+    field_names = map(x -> "field_$x", 1:NF)
+  end
+  return field_names
 end
 
 @inline function interpolate_field_values(
@@ -55,7 +91,7 @@ end
   return interps
 end
 
-# # TODO to make bcs a little clenaer
+# TODO to make bcs a little clenaer
 # @inline function map_surface_interpolants(
 #   interps::I, x_el::SVector{NxD, T}
 # ) where {I <: ReferenceFiniteElements.AbstractInterpolants, NxD, T <: Number}
@@ -64,6 +100,16 @@ end
 #   return interps
 # end
 
+# default
+function property_names(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS}
+  if NP == 0
+    prop_names = String[]
+  else
+    prop_names = map(x -> "property_$x", 1:NP)
+  end
+  return prop_names
+end
+
 """
 $(TYPEDSIGNATURES)
 """
@@ -86,6 +132,16 @@ $(TYPEDSIGNATURES)
   return SMatrix{NDof, N, T, NxNDof}(u_el)
 end
 
+# default
+function state_variable_names(::AbstractPhysics{NF, NP, NS}) where {NF, NP, NS}
+  if NS == 0
+    state_var_names = String[]
+  else
+    state_var_names = map(x -> "state_variable_$x", 1:NS)
+  end
+  return state_var_names
+end
+
 """
 $(TYPEDSIGNATURES)
 Unpacks a single fields values from a SMatrix.
@@ -121,50 +177,3 @@ crucial for performance with ```StaticArrays```.
   dofs = SVector{D, Int}(dof_start:dof_end)
   return SMatrix{D, N, T, D * N}(field[dofs, ids])
 end
-
-# Can we make something that makes interfacing with kernels easier?
-# How can we make something like this work nicely with AD?
-# struct PhysicsQuadratureState{T, ND, NN, NF, NP, NS, NDxNF, NNxND, NNxNNxND}
-#   # u::SVector{NF, T}
-#   # ∇u::SMatrix{NF, ND, T, NDxNF}
-#   # props::SVector{NP, T}
-#   # state_old::SVector{NS, T}
-#   # interpolants and gauss weight at quadrature point
-#   N::SVector{NN, T}
-#   ∇N_ξ::SMatrix{NN, ND, T, NNxND}
-#   ∇∇N_ξ::SArray{Tuple{NN, ND, ND}, T, 3, NNxNNxND}
-#   JxW::T
-#   # element level fields
-#   u_el::SMatrix{}
-# end 
-
-# physics like methods
-function damping end
-function energy end
-function mass end
-function mass! end
-function mass_action end
-function mass_action! end
-function residual end
-function residual! end
-function stiffness end
-function stiffness! end
-function stiffness_action end
-function stiffness_action! end
-
-# optimization like methods
-function gradient end
-function hessian end
-function value end
-
-# function energy(::AbstractPhysics)
-
-# end
-
-# function residual(
-#   physics, interps::ReferenceFiniteElements.Interpolants,
-#   u_el, x_el, state_old_q, props_el, dt
-# )
-#   mapped_interps = MappedInterpolants(interps, x_el)
-#   return residual(physics, mapped_interps, u_el, state_old_q, props_el, dt)
-# end
diff --git a/src/Properties.jl b/src/Properties.jl
new file mode 100644
index 0000000..38acfa8
--- /dev/null
+++ b/src/Properties.jl
@@ -0,0 +1,131 @@
+abstract type AbstractBlockProperties{V} end
+
+struct ConstantBlockProperties{
+    V <: AbstractArray{<:Number, 1}
+} <: AbstractBlockProperties{V}
+    props::V
+end
+
+struct VariableBlockProperties{
+    V <: AbstractArray{<:Number, 1}
+} <: AbstractBlockProperties{V}
+    props::V
+end
+
+struct Properties{
+    T <: Number,
+    D <: AbstractVector{T}
+} <: AbstractDiscontinuousField{T, D}
+    data::D
+    constant_block::Vector{Bool}
+    nprops::Vector{Int}
+    nepes::Vector{Int}
+    nelems::Vector{Int}
+    offsets::Vector{Int}
+    prop_names::Vector{Vector{String}}
+end
+
+# really basic case constructor where everything is constant
+# and there's a default method for all physics...
+# this is probably fragile though
+function Properties(physics)
+    data = Float64[] # TODO eventually type
+    constant_block = Bool[]
+    nprops = Int[]
+    nepes = Int[]
+    nelems = Int[]
+    offsets = Int[]
+    offset = 1
+    prop_names = Vector{String}[]
+    for physics in physics
+        props = create_properties(physics)
+        append!(data, props)
+        push!(constant_block, true)
+        push!(nprops, length(props))
+        # for this case we're just going to make em 1
+        push!(nepes, 1)
+        push!(nelems, 1)
+        push!(offsets, offset)
+        push!(prop_names, map(x -> "property_$x", 1:length(props)))
+        offset += length(props)
+    end
+    return Properties(data, constant_block, nprops, nepes, nelems, offsets, prop_names)
+end
+
+# case where everything is variable and we have a default constructor
+function Properties(fspace::FunctionSpace, physics)
+    @assert length(physics)
+    data = Float64[]
+    constant_block = Bool[]
+    nprops = Int[]
+    nepes = Int[]
+    nelems = Int[]
+    offsets = Int[]
+    offset = 1
+    prop_names = Vector{String}[]
+    @assert false "Finish me"
+end
+
+function Adapt.adapt_structure(to, props::Properties{T, D}) where {T, D}
+    data = adapt(to, props.data)
+    return Properties{T, typeof(data)}(
+        data,
+        field.constant_block,
+        field.nprops,
+        field.nepes,
+        field.nelems,
+        field.offsets
+    )
+end
+
+Base.getindex(field::Properties, i::Int) = field.data[i]
+Base.IndexStyle(::Type{<:Properties}) = IndexLinear()
+Base.size(field::Properties) = size(field.data)
+
+function Base.show(io::IO, props::Properties)
+    println(io, "Properties:")
+    for b in 1:num_blocks(props)
+        nf, nepe, ne = block_size(props, b)
+        println(io, "  Block $b:")
+        if props.constant_block[b]
+            block_props = block_view(props, b)
+            for (name, prop) in zip(props.prop_names[b], block_props.props)
+                println("    $name => $prop")
+            end
+        else
+            println(io, "    Number of properties           = $nf")
+            println(io, "    Number of entities per element = $nepe")
+            println(io, "    Number of elements             = $ne")
+            println(io, "    Property names:")
+            for name in props.prop_names
+                println(io, "      $name")
+            end
+        end
+    end
+end
+
+function Base.show(io::IO, ::MIME"text/plain", field::Properties)
+    show(io, field)
+end
+
+function block_size(field::Properties, b::Int)
+    return (field.nprops[b], field.nepes[b], field.nelems[b])
+end
+
+function block_view(props::Properties, b::Int)
+    nprops = props.nprops[b]
+    boffset = props.offsets[b]
+    if props.constant_block[b]
+        bend = boffset + nprops - 1
+        return ConstantBlockProperties(view(props.data, boffset:bend))
+    else
+        nepe = field.nepes[b]
+        nelem = field.nelems[b]
+        bend = boffset + nfield * nepe * nelem - 1
+        return VariableBlockProperties(reshape(view(field.data, boffset:bend, nprops, nepe, nelem)))
+    end
+end
+
+function num_blocks(props::Properties)
+    return length(props.nelems)
+end
diff --git a/src/Utils.jl b/src/Utils.jl
index fc5bfe6..220fe87 100644
--- a/src/Utils.jl
+++ b/src/Utils.jl
@@ -151,10 +151,11 @@ function fec_foraxes(
     end
 end
 
+# this one only does the right thing when physics is a namedtuple
 @generated function foreach_block(
     f,
     fspace::FunctionSpace{IT, IV, C, R},
-    p::AbstractParameters
+    p::Parameters
 ) where {IT, IV, C, R}
     stmts = Expr(:block)
     # NOTE need to grab one of the NamedTuple types to get the field count
@@ -164,11 +165,35 @@ end
         push!(stmts.args, quote
             f(
                 values(p.physics)[$k], values(p.properties)[$k],
-                # values(fspace.ref_fes)[$k], $k
-                fspace.ref_fes[$k], $k
-                # block_reference_element(fspace, $k), $k
+                block_reference_element(fspace, $k), $k
             )
         end)
     end
     return stmts
 end
+
+@generated function foreach_block(
+    f,
+    fspace::FunctionSpace{B, IT, IV, BTRE, C, R},
+    p::TypeStableParameters
+) where {B, IT, IV, BTRE, C, R}
+    exprs = map(1:MAX_BLOCKS) do i
+        # For each slot i, generate a branch over ALL possible ref_fe indices
+        # block_to_ref_fe_id[i] == -1 means inactive block
+        # Otherwise it's 1..length(REFS) — enumerate them all statically
+        n_refs = length(BTRE.parameters)
+        ref_dispatches = map(1:n_refs) do j
+            quote
+                if fspace.block_to_ref_fe_id[$i] == $j
+                    f(p.physics[$i], p.properties[$i], fspace.ref_fes[$j], $i)
+                end
+            end
+        end
+        quote
+            if fspace.block_to_ref_fe_id[$i] != -1
+                $(ref_dispatches...)
+            end
+        end
+    end
+    quote $(exprs...) end
+end
\ No newline at end of file
diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl
index b322923..42917ef 100644
--- a/src/assemblers/Assemblers.jl
+++ b/src/assemblers/Assemblers.jl
@@ -258,7 +258,7 @@ function _sparse_matrix(asm::AbstractAssembler, storage)
   elseif isa(type, CSRMatrix)
     return _csr_matrix(asm.matrix_pattern, asm.dof, storage)
   else
-    @assert false "Should never happen. Handle this case better"
+    @assert false "Should never happen. Handle this case better. Type is $(_sparse_matrix_type(asm))"
   end
 end
 
@@ -325,9 +325,9 @@ $(TYPEDSIGNATURES)
 assumes assemble_vector! has already been called
 """
 function residual(asm::AbstractAssembler; use_sparse_vector = false)
-  if use_sparse_vector
-    return sparsevec(asm.vector_pattern, asm.residual_unknowns)
-  else
+  # if use_sparse_vector
+    # return sparsevec(asm.vector_pattern, asm.residual_unknowns)
+  # else
     if _is_condensed(asm.dof)
       _adjust_vector_entries_for_constraints!(
         asm.residual_storage, asm.constraint_storage
@@ -341,7 +341,7 @@ function residual(asm::AbstractAssembler; use_sparse_vector = false)
       )
       return asm.residual_unknowns
     end
-  end
+  # end
 end
 
 """
diff --git a/src/bcs/DirichletBCs.jl b/src/bcs/DirichletBCs.jl
index 9136eff..8c0f078 100644
--- a/src/bcs/DirichletBCs.jl
+++ b/src/bcs/DirichletBCs.jl
@@ -347,7 +347,7 @@ function Adapt.adapt_structure(to, bcs::DirichletBCs{BCF, IV, RV}) where {BCF, I
   else
     temp_int = adapt(to, zeros(Int, 0))
     temp_floats = adapt(to, zeros(Float64, 0))
-    bc_cache = DirichletBCContainer{typeof(temp_int), typeof(temp_float)}(
+    bc_cache = DirichletBCContainer(
       temp_int, copy(temp_int),
       temp_floats, copy(temp_floats), copy(temp_floats)
     )
diff --git a/src/bcs/NeumannBCs.jl b/src/bcs/NeumannBCs.jl
index d77efc5..ab330fb 100644
--- a/src/bcs/NeumannBCs.jl
+++ b/src/bcs/NeumannBCs.jl
@@ -27,8 +27,8 @@ Internal implementation of dirichlet BCs
 """
 struct NeumannBCContainer{
   IT <: Integer,
-  IV <: AbstractArray{IT, 1},
-  RM <: AbstractArray{<:Union{<:Number, <:SVector}, 2}
+  IV <: AbstractVector{IT},
+  RM <: AbstractMatrix{<:SVector}
 } <: AbstractWeaklyEnforcedBCContainer{IT, IV, RM}
   element_conns::Connectivity{IT, IV}
   elements::IV
@@ -44,6 +44,8 @@ struct NeumannBCContainer{
   end
 end
 
+_vals_type(::Vector{NeumannBCContainer{IT, IV, RM}}) where {IT, IV, RM} = RM
+
 function Adapt.adapt_structure(to, bc::NeumannBCContainer)
   el_conns = adapt(to, bc.element_conns)
   elements = adapt(to, bc.elements)
@@ -76,8 +78,8 @@ end
 struct NeumannBCs{
   BCFuncs,
   IT      <: Integer,
-  IV      <: AbstractArray{IT, 1},
-  RM      <: AbstractArray{<:Union{<:Number, <:SVector}, 2},
+  IV      <: AbstractVector{IT},
+  RM      <: AbstractMatrix{<:SVector},
 } <: AbstractBCs{BCFuncs}
   bc_caches::Vector{NeumannBCContainer{IT, IV, RM}}
   bc_funcs::BCFuncs
@@ -96,7 +98,8 @@ end
 # handle blocks correclty 
 function NeumannBCs(mesh, dof::DofManager, neumann_bcs::Vector{NeumannBC})
   if length(neumann_bcs) == 0
-    bc_caches = NeumannBCContainer{Int, Vector{Int}, Matrix{Float64}}[]
+    ND = size(dof, 1)
+    bc_caches = NeumannBCContainer{Int, Vector{Int}, Matrix{SVector{ND, Float64}}}[]
     return NeumannBCs(bc_caches, NeumannBCFunction[], Int[], String[])
   end
 
@@ -117,8 +120,8 @@ function Adapt.adapt_structure(to, bcs::NeumannBCs)
     bc_caches = map(x -> adapt(to, x), bcs.bc_caches)
   else
     temp_int = adapt(to, zeros(Int, 0))
-    temp_floats = adapt(to, zeros(Float64, 0, 0))
-    bc_caches = NeumannBCContainer{Int, typeof(temp_int), typeof(temp_floats)}[]
+    temp_vals = adapt(to, _vals_type(bcs.bc_caches)(undef, 0, 0))
+    bc_caches = NeumannBCContainer{Int, typeof(temp_int), typeof(temp_vals)}[]
   end
   return NeumannBCs(
     bc_caches,
diff --git a/src/bcs/RobinBCs.jl b/src/bcs/RobinBCs.jl
index 92f1f2a..865bf54 100644
--- a/src/bcs/RobinBCs.jl
+++ b/src/bcs/RobinBCs.jl
@@ -27,9 +27,9 @@ Internal implementation of dirichlet BCs
 """
 struct RobinBCContainer{
   IT  <: Integer,
-  IV  <: AbstractArray{IT, 1},
-  RM  <: AbstractArray{<:Union{<:Number, <:SVector}, 2},
-  dRM <: AbstractArray{<:Union{<:Number, <:SMatrix}, 2}
+  IV  <: AbstractVector{IT},
+  RM  <: AbstractMatrix{<:SVector},
+  dRM <: AbstractMatrix{<:SMatrix}
 } <: AbstractWeaklyEnforcedBCContainer{IT, IV, RM}
   element_conns::Connectivity{IT, IV}
   elements::IV
@@ -46,6 +46,10 @@ struct RobinBCContainer{
   end
 end
 
+_dvals_type(::Vector{RobinBCContainer{IT, IV, RM, dRM}}) where {IT, IV, RM, dRM} = dRM
+_vals_type(::Vector{RobinBCContainer{IT, IV, RM, dRM}}) where {IT, IV, RM, dRM} = RM
+
+
 function Adapt.adapt_structure(to, bc::RobinBCContainer)
   el_conns = adapt(to, bc.element_conns)
   elements = adapt(to, bc.elements)
@@ -84,9 +88,9 @@ end
 struct RobinBCs{
     BCFuncs,
     IT       <: Integer,
-    IV       <: AbstractArray{IT, 1},
-    RM       <: AbstractArray{<:Union{<:Number, <:SVector}, 2},
-    dRM      <: AbstractArray{<:Union{<:Number, <:SMatrix}, 2}  
+    IV       <: AbstractVector{IT},
+    RM       <: AbstractMatrix{<:SVector},
+    dRM      <: AbstractMatrix{<:SMatrix}  
 } <: AbstractBCs{BCFuncs}
     bc_caches::Vector{RobinBCContainer{IT, IV, RM, dRM}}
     bc_funcs::BCFuncs
@@ -96,7 +100,8 @@ end
 
 function RobinBCs(mesh, dof::DofManager, robin_bcs::Vector{RobinBC})
     if length(robin_bcs) == 0
-      bc_caches = RobinBCContainer{Int, Vector{Int}, Matrix{Float64}, Matrix{Float64}}[]
+      ND = size(dof, 1)
+      bc_caches = RobinBCContainer{Int, Vector{Int}, Matrix{SVector{ND, Float64}}, Matrix{SMatrix{ND, ND, Float64, ND * ND}}}[]
       return RobinBCs(bc_caches, RobinBCFunction[], Int[], String[])
     end
 
@@ -116,8 +121,9 @@ function Adapt.adapt_structure(to, bcs::RobinBCs)
     bc_caches = map(x -> adapt(to, x), bcs.bc_caches)
   else
     temp_int = adapt(to, zeros(Int, 0))
-    temp_floats = adapt(to, zeros(Float64, 0, 0))
-    bc_caches = RobinBCContainer{Int, typeof(temp_int), typeof(temp_floats), typeof(temp_floats)}[]
+    temp_vals = adapt(to, _vals_type(bcs.bc_caches)(undef, 0, 0))
+    temp_dvals = adapt(to, _dvals_type(bcs.bc_caches)(undef, 0, 0))
+    bc_caches = RobinBCContainer{Int, typeof(temp_int), typeof(temp_vals), typeof(temp_dvals)}[]
   end
   return RobinBCs(
     bc_caches,
diff --git a/src/bcs/Sources.jl b/src/bcs/Sources.jl
index edda417..9916676 100644
--- a/src/bcs/Sources.jl
+++ b/src/bcs/Sources.jl
@@ -33,7 +33,7 @@ volume quadrature points for one block.
   and have assemblers just ping the correct block
 """
 struct SourceContainer{
-  RM <: AbstractArray{<:Union{<:Number, <:SVector}, 2}
+  RM <: AbstractMatrix{<:SVector}
 }
   vals::RM                    # size (NQ_cell, nelem) of SVector{ND}
   is_constant::Bool           # skip re-evaluation after first call
@@ -83,7 +83,7 @@ Collection of body force containers, one per specified body force entry.
 """
 struct Sources{
   SourceFuncs,
-  RM           <: AbstractArray{<:Union{<:Number, <:SVector}, 2},
+  RM           <: AbstractMatrix{<:SVector},
 } <: AbstractBCs{SourceFuncs}
   source_block_ids::Vector{Int} # note this is the id order in FEC not the id in exodus
   source_block_names::Vector{String}
@@ -102,7 +102,7 @@ struct Sources{
     caches = SourceContainer{Matrix{SVector{ND, RT}}}[]
     funcs = SourceFunction[]
     if length(sources) == 0
-      return Sources{typeof(funcs), Matrix{SVector{ND, RT}}}(Int[], String[], SourceContainer{Matrix{RT}}[], funcs)
+      return Sources{typeof(funcs), Matrix{SVector{ND, RT}}}(Int[], String[], SourceContainer{Matrix{SVector{ND, RT}}}[], funcs)
     end
 
     fspace = function_space(dof)
@@ -138,11 +138,11 @@ struct Sources{
 end
 
 function Adapt.adapt(to, sources::Sources{SourceFuncs, RM}) where {SourceFuncs, RM}
-  # needed due to failures on 1.10 and 1.11
-  if length(sources.source_caches) == 0
-    caches = Vector{SourceContainer{Matrix{Float64}}}(undef, 0)
-  else
+  if length(sources.source_caches) > 0
     caches = map(x -> adapt(to, x), sources.source_caches)
+  else
+    temp = adapt(to, _vals_type(sources.source_caches)(undef, 0, 0))
+    caches = SourceContainer{typeof(temp)}[]
   end
 
   return Sources{SourceFuncs, _vals_type(caches)}(
diff --git a/test/TestAppTools.jl b/test/TestAppTools.jl
index 8c0287c..6ad483b 100644
--- a/test/TestAppTools.jl
+++ b/test/TestAppTools.jl
@@ -132,7 +132,7 @@ end
     @test isdir("MyApp/src")
     @test isfile("MyApp/src/MyApp.jl")
     # TODO re-enable after future release
-    # AT.build_app(; path = "MyApp")
+    AT.build_app(; path = "MyApp")
     rm("MyApp"; force = true, recursive = true)
 
     mkdir("app_dir")
diff --git a/test/TestAssemblers.jl b/test/TestAssemblers.jl
index 3c63a7e..6636218 100644
--- a/test/TestAssemblers.jl
+++ b/test/TestAssemblers.jl
@@ -64,7 +64,7 @@ end
   mesh_file = "./poisson/multi_block_mesh_quad4_tri3.g"
   mesh = UnstructuredMesh(mesh_file)
   V = FunctionSpace(mesh, H1Field, Lagrange)
-  physics = Mechanics(1.0, PlaneStrain())
+  physics = Mechanics(PlaneStrain())
   props = create_properties(physics)
   u = VectorFunction(V, "displ")
   fixed(_, _) = 0.
diff --git a/test/TestMesh.jl b/test/TestMesh.jl
index 51d8046..5ec4ec8 100644
--- a/test/TestMesh.jl
+++ b/test/TestMesh.jl
@@ -1,9 +1,3 @@
-using FiniteElementContainers
-using Test
-
-geo_file_tri3 = Base.source_dir() * "/gmsh/square_meshed_with_tris.geo"
-msh_file_tri3 = Base.source_dir() * "/gmsh/square_meshed_with_tris.msh"
-
 # TODO revive this guy
 function test_amr_mesh()
   mesh = StructuredMesh("tri", (0., 0.), (1., 1.), (5, 5))
@@ -68,7 +62,7 @@ function test_amr_mesh()
   FiniteElementContainers.write_to_file(amr, "atri.exo"; force = true)
 end
 
-@testitem "Meshes - test_gmsh_mesh" begin
+@testitem "Meshes - test_gmsh_mesh" tags=[:gmsh] begin
   using Gmsh
   geo_file_tri3 = Base.source_dir() * "/gmsh/square_meshed_with_tris.geo"
   msh_file_tri3 = Base.source_dir() * "/gmsh/square_meshed_with_tris.msh"
diff --git a/test/TestProperties.jl b/test/TestProperties.jl
new file mode 100644
index 0000000..4e9a403
--- /dev/null
+++ b/test/TestProperties.jl
@@ -0,0 +1,25 @@
+@testitem "Properties - test_properties" begin
+    import FiniteElementContainers: block_size
+    using StaticArrays
+    include("mechanics/TestMechanicsCommon.jl")
+    include("poisson/TestPoissonCommon.jl")
+  
+    f(_, _) = 0.0
+    physics = (; 
+      block_1 = Poisson(f),
+      block_2 = Poisson(f)
+    )
+    props = Properties(physics)
+    display(props)
+    @test block_size(props, 1) == (0, 1, 1)
+    @test block_size(props, 2) == (0, 1, 1)
+  
+    physics = (;
+      block_1 = Mechanics(PlaneStrain()),
+      block_2 = Mechanics(PlaneStrain())
+    )
+    props = Properties(physics)
+    display(props)
+    @test block_size(props, 1) == (3, 1, 1)
+    @test block_size(props, 2) == (3, 1, 1)
+  end
\ No newline at end of file
diff --git a/test/input-file.toml b/test/input-file.toml
index 1a23e2d..1de53ac 100644
--- a/test/input-file.toml
+++ b/test/input-file.toml
@@ -17,10 +17,10 @@ file_type = "exodus"
 
 [[initial_conditions]]
 function = "one"
-block = "block_1"
-variable = "u"
+blocks = ["block_1"]
+variables = ["u"]
 
 [[boundary_conditions.dirichlet]]
 function = "zero"
-side_set = "sset_1"
-variable = "u"
+side_sets = ["sset_1"]
+variables = ["u"]
diff --git a/test/mechanics/TestMechanics.jl b/test/mechanics/TestMechanics.jl
index 1c06ef1..9d1b08f 100644
--- a/test/mechanics/TestMechanics.jl
+++ b/test/mechanics/TestMechanics.jl
@@ -21,7 +21,7 @@
 # )
   mesh = UnstructuredMesh(mesh_file)
   V = FunctionSpace(mesh, H1Field, Lagrange) 
-  physics = Mechanics(1.0, PlaneStrain())
+  physics = Mechanics(PlaneStrain())
   props = create_properties(physics)
   u = VectorFunction(V, "displ")
   times = TimeStepper(0., 1., 1)
diff --git a/test/mechanics/TestMechanicsCommon.jl b/test/mechanics/TestMechanicsCommon.jl
index a6df2fd..7b6811c 100644
--- a/test/mechanics/TestMechanicsCommon.jl
+++ b/test/mechanics/TestMechanicsCommon.jl
@@ -1,20 +1,20 @@
 # using Enzyme
 
-struct Mechanics{Form} <: AbstractPhysics{2, 2, 0}
-  density::Float64
+struct Mechanics{Form} <: AbstractPhysics{2, 3, 0}
   formulation::Form
 end
 
 function FiniteElementContainers.create_properties(::Mechanics)
+  ρ = 1e3
   K = 10.e9
   G = 1.e9
-  return SVector{2, Float64}(K, G)
+  return SVector{3, Float64}(ρ, K, G)
 end
 
 @inline function strain_energy(
   ∇u, state_old, props, dt
 )
-  K, G = props[1], props[2]
+  K, G = props[2], props[3]
   ε = symmetric(∇u)
   ψ = 0.5 * K * tr(ε)^2 + G * dcontract(dev(ε), dev(ε))
 end
@@ -22,7 +22,7 @@ end
 @inline function pk1_stress(
   ∇u, state_old, props, dt
 )
-  K, G = props[1], props[2]
+  K, G = props[2], props[3]
   ε = symmetric(∇u)
   # ε_dev = dev(ε)
   I = one(Tensor{2, 3, Float64, 9})
@@ -81,7 +81,7 @@ end
     end
   end
   M_el = SMatrix{NDOF, NDOF, eltype(N), NDOF * NDOF}(tup.data)
-  return JxW * physics.density * M_el
+  return JxW * props_el[1] * M_el
 end
 
 @inline function FiniteElementContainers.mass!(
@@ -92,7 +92,7 @@ end
 )
   interps = map_interpolants(interps, x_el)
   (; X_q, N, ∇N_X, JxW) = interps
-  scatter_with_values_and_values!(storage, physics.formulation, e, conn, N, JxW * physics.density)
+  scatter_with_values_and_values!(storage, physics.formulation, e, conn, N, JxW * props_el[1])
 end
 
 @inline function FiniteElementContainers.mass_action!(
@@ -103,7 +103,7 @@ end
 )
   interps = map_interpolants(interps, x_el)
   (; N, JxW) = interps
-  scatter_with_values_and_values!(storage, physics.formulation, e, conn, N, JxW * physics.density, v_el)
+  scatter_with_values_and_values!(storage, physics.formulation, e, conn, N, JxW * props_el[1], v_el)
 end
 
 # note for CUDA things crash without inline
diff --git a/test/poisson/TestPoisson.jl b/test/poisson/TestPoisson.jl
index 611ad49..6211b05 100644
--- a/test/poisson/TestPoisson.jl
+++ b/test/poisson/TestPoisson.jl
@@ -209,7 +209,7 @@ end
 end
 
 # MOVE below two to standalone gmsh ci testing maybe?
-@testitem "Regression test - test_poisson_dirichlet_with_nodesets_gmsh_geo_tri3" setup=[PoissonRegressionHelper] begin
+@testitem "Regression test - test_poisson_dirichlet_with_nodesets_gmsh_geo_tri3" setup=[PoissonRegressionHelper] tags=[:gmsh] begin
   using Gmsh
   geo_file_tri3 = dirname(Base.source_dir()) * "/gmsh/square_meshed_with_tris.geo"
   output_file = "poisson_test_3.e"
@@ -265,7 +265,7 @@ end
   end
 end
 
-@testitem "Regression test - test_poisson_dirichlet_with_nodesets_gmsh_msh_tri3" setup=[PoissonRegressionHelper] begin
+@testitem "Regression test - test_poisson_dirichlet_with_nodesets_gmsh_msh_tri3" setup=[PoissonRegressionHelper] tags=[:gmsh] begin
   using Gmsh
   msh_file_tri3 = dirname(Base.source_dir()) * "/gmsh/square_meshed_with_tris.msh"
   output_file = "poisson_test_4.e"
diff --git a/test/runtests.jl b/test/runtests.jl
index b7b31cd..5bc36a0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,5 +7,5 @@ using TestItems
 if "--test-amdgpu" in ARGS || "--test-cuda" in ARGS
     @run_package_tests
 else
-    @run_package_tests filter=ti->!(:gpu in ti.tags)
+    @run_package_tests filter = ti -> !(:gpu in ti.tags)
 end

From bccffedfeb0c378fd899479b44a0adb023c9b29a Mon Sep 17 00:00:00 2001
From: "Craig M. Hamel" <cmhamel32@gmail.com>
Date: Fri, 24 Apr 2026 20:25:04 -0400
Subject: [PATCH 2/2] disabling juliac testing for now.

---
 test/TestAppTools.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/TestAppTools.jl b/test/TestAppTools.jl
index 6ad483b..8c0287c 100644
--- a/test/TestAppTools.jl
+++ b/test/TestAppTools.jl
@@ -132,7 +132,7 @@ end
     @test isdir("MyApp/src")
     @test isfile("MyApp/src/MyApp.jl")
     # TODO re-enable after future release
-    AT.build_app(; path = "MyApp")
+    # AT.build_app(; path = "MyApp")
     rm("MyApp"; force = true, recursive = true)
 
     mkdir("app_dir")