diff --git a/src/Enzyme.jl b/src/Enzyme.jl index 651fbe3..e6cc8a8 100644 --- a/src/Enzyme.jl +++ b/src/Enzyme.jl @@ -56,7 +56,7 @@ function _assemble_scalar_block_enzyme_safe!( for e in axes(state_old, 3) conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) # val_el = _element_scratch(return_type, ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) @@ -168,7 +168,7 @@ function _assemble_vector_block_enzyme_safe!( conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) # # val_el = _element_scratch(return_type, ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) diff --git a/src/Fields.jl b/src/Fields.jl index 284038e..3d2841f 100644 --- a/src/Fields.jl +++ b/src/Fields.jl @@ -552,14 +552,24 @@ function properties(field::PropertyField, e::Int, b::Int) nfields = num_fields(field, b) if field.isblockconstant[b] == PROPS_CONST start = offset - finish = offset + nfields - 1 elseif field.isblockconstant[b] == PROPS_ELEMS start = offset + nfields * (e - 1) - finish = offset + nfields * e - 1 - else - @assert false "Should never happen" end - return view(field.data, start:finish) + return PropertyFieldView(field.data, start, nfields) +end + +struct PropertyFieldView{T, D <: AbstractVector{T}} <: AbstractVector{T} + data::D + start::Int + len::Int +end + +Base.size(v::PropertyFieldView) = (v.len,) +Base.length(v::PropertyFieldView) = v.len +Base.IndexStyle(::Type{<:PropertyFieldView}) = IndexLinear() +Base.@propagate_inbounds function Base.getindex(v::PropertyFieldView, i::Int) + @boundscheck checkbounds(v, i) + return @inbounds v.data[v.start + i - 1] end ###################################################################################################### diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index dfd2754..e6580b1 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -199,7 +199,11 @@ $(TYPEDSIGNATURES) offset = field.offsets[b] ec = ifelse(field.isblockconstant[b] == PROPS_CONST, 1, e) base = offset + NP * (ec - 1) - return SVector{NP, eltype(field)}(ntuple(i -> field.data[base + i - 1], NP)) + # return SVector{NP, eltype(field)}(ntuple(i -> field.data[base + i - 1], NP)) + nprops = num_fields(field, b) + # return view(field.data, base:base + NP - 1) + # return PropertyFieldView(field.data, base, NP) + return PropertyFieldView(field.data, base, nprops) end @inline function _element_level_properties(props::AbstractArray, ::Int) @@ -431,7 +435,7 @@ function _assemble_block!( fec_foraxes(state_old, 3) do e conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) val_el = _element_scratch(return_type, ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) @@ -463,7 +467,7 @@ function _assemble_block!( fec_foraxes(state_old, 3) do e conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old = element_level_fields(ref_fe, conn, X, U, U_old) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) state_old_q = _quadrature_level_state(state_old, q, e) diff --git a/src/assemblers/MatrixAction.jl b/src/assemblers/MatrixAction.jl index 5dc1ea7..e45d3cd 100644 --- a/src/assemblers/MatrixAction.jl +++ b/src/assemblers/MatrixAction.jl @@ -66,7 +66,7 @@ function _assemble_block_matrix_free_action!( fec_foraxes(state_old, 3) do e conn = connectivity(ref_fe, conns, e, coffset) x_el, u_el, u_el_old, v_el = element_level_fields(ref_fe, conn, X, U, U_old, V) - props_el = properties(props, physics, e, b) + props_el = properties(props, e, b) Kv_el = _element_scratch(AssembledVector(), ref_fe, U) for q in 1:num_cell_quadrature_points(ref_fe) interps = _cell_interpolants(ref_fe, q) diff --git a/test/TestFields.jl b/test/TestFields.jl index 96aa5a0..847d453 100644 --- a/test/TestFields.jl +++ b/test/TestFields.jl @@ -281,6 +281,41 @@ end @test all(FiniteElementContainers.properties(props, 1, 1) .≈ original) end +@testitem "Fields - test_property_field_view_eltype_is_concrete" begin + using StaticArrays + # `PropertyFieldView` must carry the element type as a parameter. Declaring + # it as `AbstractVector{eltype(D)}` over a bare `D` silently yields + # `AbstractVector{Any}`, because `eltype` of an unbound TypeVar is `Any`. + # Downstream that is not cosmetic: ConstitutiveModels' `module_props` builds + # `SVector{NP, eltype(props)}`, so an `Any` eltype turns every constitutive + # evaluation into a boxed, dynamically dispatched call. + props = FiniteElementContainers.PropertyField([[1.0, 2.0, 3.0]]) + view = FiniteElementContainers.properties(props, 1, 1) + @test eltype(view) === Float64 + @test view isa AbstractVector{Float64} + @test eltype(collect(view)) === Float64 + @test isconcretetype(eltype(SVector{2, eltype(view)}(view[1], view[2]))) + @test Base.IndexStyle(typeof(view)) === IndexLinear() +end + +@testitem "Fields - test_property_field_view_is_bounds_checked" begin + # Every block's properties live in one flat vector, so an unchecked + # out-of-range read returns the *next* block's properties instead of + # failing. Two blocks with different property counts make that concrete. + props = FiniteElementContainers.PropertyField([[1.0, 2.0, 3.0], [10.0, 20.0, 30.0, 40.0]]) + block1 = FiniteElementContainers.properties(props, 1, 1) + @test length(block1) == 3 + @test block1[3] ≈ 3.0 + # Without a bounds check this returns 10.0 -- block 2's first property. + @test_throws BoundsError block1[4] + @test_throws BoundsError block1[0] + + block2 = FiniteElementContainers.properties(props, 1, 2) + @test length(block2) == 4 + @test block2[4] ≈ 40.0 + @test_throws BoundsError block2[5] +end + @testitem "Fields - test_state_variable_field" begin a1 = rand(2, 3, 40) a2 = rand(3, 4, 10)