From fd585934f53eaf2379af1ff5983015b614ecea55 Mon Sep 17 00:00:00 2001 From: Alejandro Mota Date: Sat, 25 Jul 2026 01:56:24 -0700 Subject: [PATCH] Give element blocks one canonical order, and match physics by name Per-block mesh data lives in `Dict`s keyed by block name, while the block names themselves live in an ordered `Vector`. `Dict` iteration is hash order, so `values(mesh.element_conns)` and `mesh.element_block_names` are two different orderings of the same blocks, and code that indexes blocks positionally was mixing them: * `FunctionSpace` took connectivity and element id maps from `values(...)` but block names and reference elements from `element_block_names`. On a three block mesh named b1/b2/b3 with 1, 2 and 3 elements, `block_names(fspace)` read ["b1","b2","b3"] while `num_elements(fspace, b)` read [3, 2, 1] -- block `b`'s name and block `b`'s elements were different blocks. * `BCBookKeeping` and `_setup_sideset` numbered blocks with `enumerate(values(mesh.element_id_maps))` and then used that number to index `element_block_names`, so a side set resolved to the wrong block name, and from there to the wrong connectivity and reference element. * `write_to_file` wrote block names in file order but the blocks themselves in `Dict` order, pairing each block with another's name. Introduce `block_names(mesh)` as the single source of truth for what "block b" means, with `block_conns` / `block_id_maps` returning per-block data in that same order, and route every positional use through them. `block_names(fspace)` exposes the same order downstream. With one order established, `physics` and `properties` can be matched to blocks by name instead of by position. `Parameters` now permutes a supplied `NamedTuple` into block order and requires its keys to be exactly the block names, replacing the `TODO re-arrange physics tuple to match fspaces` -- previously a `NamedTuple` written in any other order silently gave each block another block's material, with nothing to flag it. A single physics/properties object is still shared across all blocks, but is now keyed by the real block names rather than invented `region_N` placeholders. Nothing caught any of this because every multi-block fixture in the suite has two blocks named block_1/block_2, whose hash order happens to equal their file order. `test/TestBlockOrdering.jl` adds a three block fixture where it does not, with blocks of different sizes so a permutation shows up as a size and not just a name. Behavior changes for downstream users: `p.physics` / `p.properties` are keyed by block name rather than `region_N`, and a `NamedTuple` whose keys are not exactly the block names now raises `BlockMismatchError` instead of being accepted. --- Project.toml | 2 +- docs/src/meshes.md | 27 +++++ docs/src/parameters.md | 37 +++++++ src/FiniteElementContainers.jl | 1 + src/FunctionSpaces.jl | 44 +++++--- src/Parameters.jl | 78 +++++++++++--- src/bcs/BoundaryConditions.jl | 27 +++-- src/bcs/Sources.jl | 4 +- src/meshes/AMRMesh.jl | 6 ++ src/meshes/Meshes.jl | 41 ++++++-- test/TestBlockOrdering.jl | 183 +++++++++++++++++++++++++++++++++ 11 files changed, 397 insertions(+), 53 deletions(-) create mode 100644 test/TestBlockOrdering.jl diff --git a/Project.toml b/Project.toml index c2bc4040..af01ae05 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FiniteElementContainers" uuid = "d08262e4-672f-4e7f-a976-f2cea5767631" -version = "0.14.3" +version = "0.15.0" authors = ["Craig M. Hamel and contributors"] [deps] diff --git a/docs/src/meshes.md b/docs/src/meshes.md index 0d9d99b8..e396021a 100644 --- a/docs/src/meshes.md +++ b/docs/src/meshes.md @@ -205,6 +205,33 @@ returns the connectivity and element type associated with the block named `"soli This organization makes it straightforward to assign different material models or physics objects to different regions of a mesh. +## Block order + +Per-block mesh data is stored in `Dict`s keyed by block name, but a great deal +of the library refers to blocks by *position* — `num_elements(fspace, b)`, +`foreach_block`, and the per-block entries of `physics` and `properties` all +index block `b`. `Dict` iteration is hash order, which is neither the order the +blocks appear in the mesh file nor stable under a change of block names, so it +must never be used to establish that position. + +The canonical order is given by + +```julia +block_names(mesh) +``` + +and block `b` is `block_names(mesh)[b]`. The matching per-block data in the same +order is available as + +```julia +block_conns(mesh) +block_id_maps(mesh) +``` + +A `FunctionSpace` carries the same order, so `block_names(fspace)[b]` names the +block whose connectivity, reference element, physics and properties all live at +index `b`. + --- # Nodesets and Sidesets diff --git a/docs/src/parameters.md b/docs/src/parameters.md index 672cd0d0..4176969d 100644 --- a/docs/src/parameters.md +++ b/docs/src/parameters.md @@ -1,4 +1,41 @@ # Parameters + +`Parameters` bundles everything a physics evaluation needs that is not the +solution field itself: boundary and initial conditions, sources, the time +stepper, and the per-block `physics` and `properties`. + +## Per-block physics and properties + +`physics` and `properties` may be given either as a single object shared by the +whole mesh + +```julia +p = create_parameters(mesh, asm, physics, props) +``` + +or as a `NamedTuple` with one entry per element block, keyed by block name + +```julia +p = create_parameters( + mesh, asm, + (steel = steel_physics, foam = foam_physics), + (steel = steel_props, foam = foam_props) +) +``` + +In the single-object form the object is replicated across every block. In the +`NamedTuple` form the keys must be exactly the mesh's block names: a block with +no entry, or an entry naming no block, raises a `BlockMismatchError` rather than +running with some block silently taking another block's material. + +The entries are reordered to match the block order of the function space (see +[Block order](@ref)), so the order the `NamedTuple` is written in does not +matter — only the names do. Downstream, `p.physics` and `p.properties` are +always keyed by block name and ordered by block index, which is the pairing +`foreach_block` and the assembly kernels rely on. + +## API + ```@autodocs Modules = [FiniteElementContainers] Pages = ["Parameters.jl"] diff --git a/src/FiniteElementContainers.jl b/src/FiniteElementContainers.jl index da33f4f9..8d9b541f 100644 --- a/src/FiniteElementContainers.jl +++ b/src/FiniteElementContainers.jl @@ -110,6 +110,7 @@ export evolve! export FileMesh export StructuredMesh export UnstructuredMesh +export block_names export distribute_mesh export element_blocks export element_ids diff --git a/src/FunctionSpaces.jl b/src/FunctionSpaces.jl index 4234b9ea..b7f2e3d3 100644 --- a/src/FunctionSpaces.jl +++ b/src/FunctionSpaces.jl @@ -65,25 +65,25 @@ 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) - for block_name in mesh.element_block_names + for block_name in block_names(mesh) el_type = mesh.element_types[block_name] push!(block_ids, _el_name_to_juliac_safe_id[el_type]) end return block_ids else - return 1:length(mesh.element_types) |> collect + return 1:length(block_names(mesh)) |> collect end end function _setup_juliac_safe_block_to_ref_fe_id(mesh::AbstractMesh) - names = mesh.element_block_names + names = block_names(mesh) 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 + return 1:length(block_names(mesh)) |> collect end # Lagrange elements @@ -102,13 +102,13 @@ const _juliac_safe_ref_fes = ( default code path that sets up ref fes as a namedtuple """ function _setup_ref_fes( - mesh::AbstractMesh, + mesh::AbstractMesh, interp_type, p_degree, q_type::Type{<:ReferenceFiniteElements.AbstractQuadratureType}, q_degree ) - block_names = mesh.element_block_names + names = block_names(mesh) ref_fes = ReferenceFE[] - for block_name in block_names + for block_name in names elem_name = mesh.element_types[block_name] elem_type = elem_type_map[uppercase(elem_name)] if p_degree === nothing @@ -121,7 +121,7 @@ function _setup_ref_fes( ref_fe = ReferenceFE(elem_type{interp_type, p_degree}(), q_type(q_degree)) push!(ref_fes, ref_fe) end - ref_fes = NamedTuple{tuple(Symbol.(values(block_names))...)}(tuple(ref_fes...)) + ref_fes = NamedTuple{tuple(Symbol.(names)...)}(tuple(ref_fes...)) return ref_fes end @@ -220,9 +220,11 @@ function FunctionSpace{is_juliac_safe}( ref_fes = _setup_ref_fes(mesh, interp_type, nothing, q_type, q_degree) end coords = mesh.nodal_coords - conns = Connectivity([val for val in values(mesh.element_conns)]) + # canonical block order -- NOT `values(mesh.element_conns)`, whose Dict + # iteration order does not agree with `block_names(mesh)` + conns = Connectivity(block_conns(mesh)) end - elem_id_maps = [val for val in values(mesh.element_id_maps)] + elem_id_maps = block_id_maps(mesh) if is_juliac_safe block_to_ref_fe_id = _setup_juliac_safe_block_to_ref_fe_id(mesh) else @@ -231,7 +233,7 @@ function FunctionSpace{is_juliac_safe}( end return FunctionSpace{is_juliac_safe}( - mesh.element_block_names, block_to_ref_fe_id, coords, + block_names(mesh), block_to_ref_fe_id, coords, conns, elem_id_maps, mesh.node_id_map, ref_fes ) end @@ -249,22 +251,22 @@ function FunctionSpace{is_juliac_safe}( ref_fes = _setup_ref_fes(mesh, interp_type, p_degree, q_type, q_degree) end - conns = Connectivity([val for val in values(mesh.element_conns)]) - coords = L2Field(map(x -> mesh.nodal_coords[:, x], [values(mesh.element_conns)...])) + # canonical block order throughout -- the L2 coordinates and the offsets that + # index into them have to agree with each other AND with `block_names(mesh)` + coords = L2Field(map(x -> mesh.nodal_coords[:, x], block_conns(mesh))) new_conns = Array{Int, 2}[] offset = 1 - for name in keys(mesh.element_conns) - conn = mesh.element_conns[name] + for conn in block_conns(mesh) push!(new_conns, reshape(offset:offset + length(conn) - 1, size(conn)...)) offset += size(conn, 1) * size(conn, 2) end conns = Connectivity(new_conns) - elem_id_maps = [val for val in values(mesh.element_id_maps)] + elem_id_maps = block_id_maps(mesh) block_to_ref_fe_id = _setup_block_to_ref_fe_id(mesh, is_juliac_safe) return FunctionSpace{is_juliac_safe}( - mesh.element_block_names, block_to_ref_fe_id, coords, + block_names(mesh), block_to_ref_fe_id, coords, conns, elem_id_maps, mesh.node_id_map, ref_fes ) end @@ -294,6 +296,14 @@ function _is_juliac_safe(::FunctionSpace{B, I, V, BTRE, C, R}) where {B, I, V, B return B end +""" +$(TYPEDSIGNATURES) +Names of the element blocks, in block order: `block_names(fspace)[b]` is the +name of the block whose connectivity, reference element, physics and properties +all live at index `b`. +""" +block_names(fspace::FunctionSpace) = fspace.block_names + function block_entity_size(fspace::FunctionSpace, b::Int) return (num_entities_per_element(fspace, b), num_elements(fspace, b)) end diff --git a/src/Parameters.jl b/src/Parameters.jl index d53728c4..5e79a347 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -1,3 +1,56 @@ +struct BlockMismatchError <: AbstractFECError + msg::String +end +_block_mismatch_error(msg::String) = throw(BlockMismatchError(msg)) + +function _check_block_keys(given, expected, what) + extra = filter(x -> !(x in expected), collect(given)) + absent = filter(x -> !(x in given), collect(expected)) + if !isempty(extra) || !isempty(absent) + msg = "$what must have exactly one entry per element block.\n" * + " element blocks : $(join(expected, ", "))\n" * + " $what entries : $(join(given, ", "))" + isempty(extra) || (msg *= "\n not element blocks : $(join(extra, ", "))") + isempty(absent) || (msg *= "\n blocks with no entry : $(join(absent, ", "))") + _block_mismatch_error(msg) + end + return nothing +end + +""" +Align a user-supplied `physics`/`properties` argument with the element blocks of +`fspace`, returning a `NamedTuple` keyed by block name and ordered by block +index, so that entry `b` always belongs to block `b`. + +Everything downstream -- `_setup_state_variables`, `foreach_block`, the +assembly kernels -- pairs entry `b` with block `b` positionally. A `NamedTuple` +supplied in a different order than the mesh's blocks would therefore hand each +block another block's material without any error, so it is permuted here, and +its keys are required to be exactly the block names. +""" +function _align_blocks(fspace, x::NamedTuple, what) + names = tuple(Symbol.(block_names(fspace))...) + _check_block_keys(keys(x), names, what) + return NamedTuple{names}(map(name -> getfield(x, name), names)) +end + +# A bare `Tuple` carries no block names, so there is no way to tell whether it +# is in block order or not. Rejecting is the only safe reading. +function _align_blocks(fspace, x::Tuple, what) + _block_mismatch_error( + "$what was given as an unnamed Tuple, which cannot be matched to element " * + "blocks. Supply a NamedTuple keyed by block name " * + "($(join(block_names(fspace), ", "))), or a single value to share across " * + "all blocks." + ) +end + +# a single physics/properties object shared by every block +function _align_blocks(fspace, x, what) + names = tuple(Symbol.(block_names(fspace))...) + return NamedTuple{names}(ntuple(_ -> x, length(names))) +end + function _setup_state_variables(fspace, physics) state_old = Array{Float64, 3}[] state_new = Array{Float64, 3}[] @@ -96,22 +149,8 @@ function Parameters( end # for mixed spaces we'll need to do this more carefully - if isa(physics, AbstractPhysics) - syms = map(x -> Symbol("region_$x"), 1:length(fspace.ref_fes)) - physics = map(x -> physics, syms) - physics = NamedTuple{tuple(syms...)}(tuple(physics...)) - else - @assert isa(physics, NamedTuple) - # TODO re-arrange physics tuple to match fspaces when appropriate - end - - if isa(properties, AbstractArray) - syms = map(x -> Symbol("region_$x"), 1:length(fspace.ref_fes)) - properties = map(x -> properties, syms) - properties = NamedTuple{tuple(syms...)}(tuple(properties...)) - else - @assert isa(properties, NamedTuple) - end + physics = _align_blocks(fspace, physics, "physics") + properties = _align_blocks(fspace, properties, "properties") # setup state variables state_old, state_new = _setup_state_variables(fspace, physics) @@ -234,6 +273,9 @@ struct TypeStableParameters{ pbcs = PeriodicBCs{SF}(mesh, dof, pbcs) srcs = Sources{VF}(mesh, dof, srcs) + physics = _align_blocks(fspace, physics, "physics") + props = _align_blocks(fspace, props, "properties") + state_old, state_new = _setup_state_variables(fspace, physics) coords = mesh.nodal_coords @@ -260,12 +302,16 @@ struct TypeStableParameters{ ) where {SF, VF} dof = assembler.dof ND = size(dof, 1) + fspace = function_space(dof) ics = InitialConditions{SF}(mesh, dof, ics) dbcs = DirichletBCs{SF}(mesh, dof, dbcs) nbcs = NeumannBCs{VF}(mesh, dof, nbcs) pbcs = PeriodicBCs{SF}(mesh, dof, pbcs) srcs = Sources{VF}(mesh, dof, srcs) + physics = _align_blocks(fspace, physics, "physics") + props = _align_blocks(fspace, props, "properties") + coords = mesh.nodal_coords field = create_field(assembler) field_old = create_field(assembler) diff --git a/src/bcs/BoundaryConditions.jl b/src/bcs/BoundaryConditions.jl index c9d9726c..871b13d4 100644 --- a/src/bcs/BoundaryConditions.jl +++ b/src/bcs/BoundaryConditions.jl @@ -141,7 +141,11 @@ function BCBookKeeping( # gather the blocks that are present in this sideset # and also map global element id to local element id # TODO this isn't quite right - for (n, val) in enumerate(values(mesh.element_id_maps)) + # NOTE `n` is used below to index into `block_names(mesh)`, so it has to be + # a canonical block index -- iterating `values(mesh.element_id_maps)` would + # number the blocks in Dict hash order instead. + for (n, name) in enumerate(block_names(mesh)) + val = mesh.element_id_maps[name] # note these are the local elem id to the block, e.g. starting from 1. indices_in_sset = indexin(val, elements) filter!(x -> x !== nothing, indices_in_sset) @@ -152,7 +156,7 @@ function BCBookKeeping( end @assert length(unique(blocks)) == 1 "Sidesets need to be in a single block" - block_name = mesh.element_block_names[blocks[1]] + block_name = block_names(mesh)[blocks[1]] indices_in_sset = indexin(elements, mesh.element_id_maps[block_name]) filter!(x -> x !== nothing, indices_in_sset) elements = convert(Vector{Int}, indices_in_sset) @@ -274,7 +278,7 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) new_bcs = type[] new_funcs = Function[] block_ids = Int[] - block_names = String[] + bc_block_names = String[] # sideset_ids = Int[] # sideset_names = String[] for (bk, func) in zip(bks, funcs) @@ -287,9 +291,9 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) for block in blocks block_name = mesh.element_block_names_map[block] - block_id = findfirst(x -> x == block_name, mesh.element_block_names) + block_id = findfirst(x -> x == block_name, block_names(mesh)) push!(block_ids, block_id) - push!(block_names, block_name) + push!(bc_block_names, block_name) ids = findall(x -> x == block, bk.blocks) new_blocks = bk.blocks[ids] @@ -299,7 +303,7 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) # TODO update nodes and dofs new_bk = BCBookKeeping(new_blocks, bk.dofs, new_elements, bk.nodes, new_sides, new_side_nodes) - id = findfirst(x -> x == block_name, mesh.element_block_names) + id = findfirst(x -> x == block_name, block_names(mesh)) ref_fe = fspace.ref_fes[id] NQ = num_surface_quadrature_points(ref_fe) ND = length(dof.var) @@ -331,7 +335,7 @@ function _setup_weakly_enforced_bc_container(mesh, dof, bcs, type) push!(new_funcs, func) end end - return new_bcs, new_funcs, block_ids, block_names + return new_bcs, new_funcs, block_ids, bc_block_names end function _setup_sideset(mesh, dof, bc) @@ -348,7 +352,10 @@ function _setup_sideset(mesh, dof, bc) # gather the blocks that are present in this sideset # and also map global element id to local element id blocks = Vector{Int64}(undef, 0) - for (n, val) in enumerate(values(mesh.element_id_maps)) + # NOTE `n` indexes `block_names(mesh)` below, so it has to be a canonical + # block index -- see the same loop in `BCBookKeeping`. + for (n, name) in enumerate(block_names(mesh)) + val = mesh.element_id_maps[name] # note these are the local elem id to the block, e.g. starting from 1. indices_in_sset = indexin(val, elements) filter!(x -> x !== nothing, indices_in_sset) @@ -360,7 +367,7 @@ function _setup_sideset(mesh, dof, bc) unique_block_ids = sort(unique(blocks)) @assert length(unique_block_ids) == 1 "Sidesets need to be in a single block" - block_name = mesh.element_block_names[unique_block_ids[1]] + block_name = block_names(mesh)[unique_block_ids[1]] ids = findall(x -> x == unique_block_ids[1], blocks) indices_in_sset = indexin(elements, mesh.element_id_maps[block_name]) filter!(x -> x !== nothing, indices_in_sset) @@ -368,7 +375,7 @@ function _setup_sideset(mesh, dof, bc) # end bc bookkeeping - block_id = findfirst(x -> x == block_name, mesh.element_block_names) + block_id = findfirst(x -> x == block_name, block_names(mesh)) # this probably doesn't do anything for when we just have one block blocks = blocks[ids] diff --git a/src/bcs/Sources.jl b/src/bcs/Sources.jl index c5370371..66089e92 100644 --- a/src/bcs/Sources.jl +++ b/src/bcs/Sources.jl @@ -113,7 +113,7 @@ struct Sources{ # source_block_names = String[] # for source in sources # block_name = source.block_name - # block_id = findfirst(x -> x == block_name, mesh.element_block_names) + # block_id = findfirst(x -> x == block_name, block_names(mesh)) # push!(source_block_ids, block_id) # push!(source_block_names, block_name) # NQ, NE = block_quadrature_size(fspace, block_id) @@ -175,7 +175,7 @@ struct Sources{ source_block_names = String[] for source in sources block_name = source.block_name - block_id = findfirst(x -> x == block_name, mesh.element_block_names) + block_id = findfirst(x -> x == block_name, block_names(mesh)) push!(source_block_ids, block_id) push!(source_block_names, block_name) NQ, NE = block_quadrature_size(fspace, block_id) diff --git a/src/meshes/AMRMesh.jl b/src/meshes/AMRMesh.jl index d2145777..55335773 100644 --- a/src/meshes/AMRMesh.jl +++ b/src/meshes/AMRMesh.jl @@ -51,6 +51,12 @@ function AMRMesh(mesh::AbstractMesh) ) end +# `AMRMesh` stores block names as an id => name `Dict` rather than an ordered +# `Vector`, so canonical order is by ascending block id -- which is the order the +# blocks appear in the mesh file, matching every other mesh type. +block_names(mesh::AMRMesh) = + [mesh.element_block_names[id] for id in sort!(collect(keys(mesh.element_block_names)))] + # this uses a longest edge initialization # slightly smarter than simply initalizing as 1 # but not as great as structured based diff --git a/src/meshes/Meshes.jl b/src/meshes/Meshes.jl index d8d2601e..f7c6f94e 100644 --- a/src/meshes/Meshes.jl +++ b/src/meshes/Meshes.jl @@ -53,16 +53,42 @@ $(TYPEDEF) """ abstract type AbstractMesh end +""" +$(TYPEDSIGNATURES) +Element block names in canonical block order, i.e. the order the blocks appear +in the mesh file. + +This is the single source of truth for what "block `b`" means. Meshes store +per-block data (`element_conns`, `element_id_maps`, `element_types`) in `Dict`s +keyed by block name, and `Dict` iteration order is hash order -- neither the +mesh file order nor stable under a change of block names. Anything that indexes +blocks by position must derive that order from here rather than from +`values(some_dict)`. +""" +block_names(mesh::AbstractMesh) = mesh.element_block_names + +""" +$(TYPEDSIGNATURES) +Per-block element connectivities in canonical block order. See [`block_names`](@ref). +""" +block_conns(mesh::AbstractMesh) = [mesh.element_conns[name] for name in block_names(mesh)] + +""" +$(TYPEDSIGNATURES) +Per-block element ID maps in canonical block order. See [`block_names`](@ref). +""" +block_id_maps(mesh::AbstractMesh) = [mesh.element_id_maps[name] for name in block_names(mesh)] + function Base.show(io::IO, mesh::AbstractMesh) println(io, typeof(mesh).name.name, ":") println(io, " Number of dimensions = $(size(mesh.nodal_coords, 1))") println(io, " Number of nodes = $(size(mesh.nodal_coords, 2))") println(io, " Element Blocks:") - for (name, type) in zip(values(mesh.element_block_names), mesh.element_types) + for name in block_names(mesh) conn = mesh.element_conns[name] println(io, " $name:") - println(io, " Element type = $type") + println(io, " Element type = $(mesh.element_types[name])") println(io, " Number of elements = $(size(conn, 2))") end @@ -270,13 +296,14 @@ function write_to_file(mesh::AbstractMesh, file_name::String; force::Bool = fals # write_id_map(exo, NodeMap, convert.(Int32, mesh.node_id_map)) # write block names - # block_names = map(String, values(mesh.element_block_names)) - block_names = mesh.element_block_names - write_names(exo, Block, block_names) + write_names(exo, Block, block_names(mesh)) # TODO write block id maps - # for (n, block_name) in mesh.element_block_names - for (n, block_name) in mesh.element_block_names_map + # Blocks are written in ascending block id, matching the order the names were + # just written in -- iterating the id => name `Dict` directly would write them + # in hash order and pair each block with another block's name. + for n in sort!(collect(keys(mesh.element_block_names_map))) + block_name = mesh.element_block_names_map[n] el_type = mesh.element_types[block_name] conn = mesh.element_conns[block_name] write_block(exo, n, String(el_type), conn |> collect) diff --git a/test/TestBlockOrdering.jl b/test/TestBlockOrdering.jl new file mode 100644 index 00000000..e1be4477 --- /dev/null +++ b/test/TestBlockOrdering.jl @@ -0,0 +1,183 @@ +# Regression tests for element block ORDER. +# +# Meshes keep per-block data (`element_conns`, `element_id_maps`, ...) in `Dict`s +# keyed by block name, but the block names themselves live in an ordered +# `Vector`. `Dict` iteration is hash order, so `values(mesh.element_conns)` and +# `mesh.element_block_names` are two *different* orderings of the same blocks. +# Anything that indexes blocks positionally has to pick one, and the two were +# mixed: +# +# * `FunctionSpace` took connectivity from `values(mesh.element_conns)` but +# block names and reference elements from `mesh.element_block_names`, so +# `block_names(fspace)[b]` did not name the block whose elements sit at +# index `b`. +# * `BCBookKeeping` numbered blocks by `enumerate(values(mesh.element_id_maps))` +# and then used that number to index `mesh.element_block_names`, so a side +# set resolved to the wrong block name. +# * `Parameters` paired `physics`/`properties` entry `b` with block `b` without +# checking or reordering, so a `NamedTuple` written in any other order handed +# each block another block's material. +# +# None of it was caught because every multi-block fixture in the suite has two +# blocks named `block_1`/`block_2`, for which the hash order happens to equal +# the file order. Three blocks named `b1`/`b2`/`b3` do not have that property, +# which is what this fixture is for. +# +# The blocks are given DIFFERENT element counts (1, 2, 3) so that a permutation +# is visible as a size, not just as a name. + +@testsnippet BlockOrderingHelper begin + using Exodus + using StaticArrays + include("poisson/TestPoissonCommon.jl") + + # A 6 x 1 strip of QUAD4s over [0,6] x [0,1], split into three blocks holding + # 1, 2 and 3 elements respectively. + const BLOCK_NAMES = ["b1", "b2", "b3"] + const BLOCK_SIZES = [1, 2, 3] + + function write_three_block_mesh(path) + nx, ny = 6, 1 + coords = zeros(Float64, 2, (nx + 1) * (ny + 1)) + nid(i, j) = i + (j - 1) * (nx + 1) + for j in 1:(ny + 1), i in 1:(nx + 1) + coords[1, nid(i, j)] = float(i - 1) + coords[2, nid(i, j)] = float(j - 1) + end + quad(i) = Int32[nid(i, 1), nid(i + 1, 1), nid(i + 1, 2), nid(i, 2)] + + isfile(path) && rm(path) + init = Initialization( + Int32(2), Int32(size(coords, 2)), Int32(nx), + Int32(length(BLOCK_NAMES)), Int32(0), Int32(0) + ) + exo = ExodusDatabase{Int32, Int32, Int32, Float64}(path, "w", init) + write_coordinates(exo, coords) + + first_elem = 1 + for (b, (name, n)) in enumerate(zip(BLOCK_NAMES, BLOCK_SIZES)) + conn = reduce(hcat, [quad(i) for i in first_elem:(first_elem + n - 1)]) + write_block(exo, Block(Int32(b), size(conn, 2), 4, "QUAD4", conn)) + write_name(exo, Block(exo, b), name) + first_elem += n + end + close(exo) + return path + end +end + +@testitem "Block ordering - mesh and function space" setup=[BlockOrderingHelper] begin + import FiniteElementContainers: block_names, block_conns, block_id_maps + + mktempdir() do dir + mesh = UnstructuredMesh(write_three_block_mesh(joinpath(dir, "three_block.g"))) + + # The premise of this fixture. If a future Julia gives `Dict` insertion + # order, the test below still passes but no longer exercises the bug, and + # the block names should be re-chosen. + if collect(keys(mesh.element_conns)) == mesh.element_block_names + @warn "Dict order now matches file order; this fixture no longer " * + "exercises the block ordering bug -- pick different block names." + end + + @test block_names(mesh) == BLOCK_NAMES + @test [size(c, 2) for c in block_conns(mesh)] == BLOCK_SIZES + @test [length(m) for m in block_id_maps(mesh)] == BLOCK_SIZES + + V = FunctionSpace(mesh, H1Field, Lagrange) + + # The invariant: name, connectivity and reference element at index `b` all + # belong to the same block. Before the fix `block_names(V)` read + # ["b1","b2","b3"] while `num_elements(V, b)` read [3, 2, 1]. + @test block_names(V) == BLOCK_NAMES + @test [num_elements(V, b) for b in 1:length(BLOCK_NAMES)] == BLOCK_SIZES + @test collect(keys(V.ref_fes)) == Symbol.(BLOCK_NAMES) + for (b, name) in enumerate(BLOCK_NAMES) + @test num_elements(V, b) == size(mesh.element_conns[name], 2) + end + end +end + +@testitem "Block ordering - physics and properties alignment" setup=[BlockOrderingHelper] begin + import FiniteElementContainers: block_names, _align_blocks, BlockMismatchError + + mktempdir() do dir + mesh = UnstructuredMesh(write_three_block_mesh(joinpath(dir, "three_block.g"))) + V = FunctionSpace(mesh, H1Field, Lagrange) + + # A NamedTuple in any order is permuted into block order, and keeps its + # association with the block it was named for. + for order in ((b1 = 1, b2 = 2, b3 = 3), + (b3 = 3, b1 = 1, b2 = 2), + (b2 = 2, b3 = 3, b1 = 1)) + aligned = _align_blocks(V, order, "physics") + @test keys(aligned) == (:b1, :b2, :b3) + @test values(aligned) == (1, 2, 3) + end + + # A single object is shared by every block, and is keyed by the real block + # names rather than invented `region_N` placeholders. + shared = _align_blocks(V, :one_material, "physics") + @test keys(shared) == (:b1, :b2, :b3) + @test all(v -> v === :one_material, values(shared)) + + # A block with no entry, an entry naming no block, and an unnamed Tuple all + # have to be rejected: each of them would otherwise silently give some block + # another block's material. + @test_throws BlockMismatchError _align_blocks(V, (b1 = 1, b2 = 2), "physics") + @test_throws BlockMismatchError _align_blocks(V, (b1 = 1, b2 = 2, b3 = 3, b4 = 4), "physics") + @test_throws BlockMismatchError _align_blocks(V, (b1 = 1, b2 = 2, typo = 3), "physics") + @test_throws BlockMismatchError _align_blocks(V, (1, 2, 3), "physics") + + # The error has to name the blocks involved, otherwise it is no better than + # the silent misassignment it replaces. + err = try + _align_blocks(V, (b1 = 1, b2 = 2, typo = 3), "physics") + catch e + sprint(showerror, e) + end + @test occursin("typo", err) + @test occursin("b3", err) + end +end + +@testitem "Block ordering - parameters" setup=[BlockOrderingHelper] begin + import FiniteElementContainers: block_names, block_size, BlockMismatchError + + mktempdir() do dir + mesh = UnstructuredMesh(write_three_block_mesh(joinpath(dir, "three_block.g"))) + V = FunctionSpace(mesh, H1Field, Lagrange) + u = ScalarFunction(V, "u") + asm = SparseMatrixAssembler(u) + + f(X, _) = 0.0 + one_physics = Poisson(f) + one_props = create_properties(one_physics) + + # Single material for the whole mesh: replicated, and keyed by block name. + p = create_parameters(mesh, asm, one_physics, one_props) + @test collect(keys(p.physics)) == Symbol.(BLOCK_NAMES) + @test collect(keys(p.properties)) == Symbol.(BLOCK_NAMES) + + # State variables are allocated per block by walking `values(physics)` + # against `block_quadrature_size(fspace, b)`, so their element counts are a + # direct check that entry `b` really is block `b`. + @test [block_size(p.state_old, b)[3] for b in 1:length(BLOCK_NAMES)] == BLOCK_SIZES + + # Per-block materials supplied out of order still land on the right blocks. + scrambled_physics = (b3 = Poisson(f), b1 = Poisson(f), b2 = Poisson(f)) + scrambled_props = (b3 = one_props, b1 = one_props, b2 = one_props) + p = create_parameters(mesh, asm, scrambled_physics, scrambled_props) + @test collect(keys(p.physics)) == Symbol.(BLOCK_NAMES) + @test collect(keys(p.properties)) == Symbol.(BLOCK_NAMES) + @test [block_size(p.state_old, b)[3] for b in 1:length(BLOCK_NAMES)] == BLOCK_SIZES + + # And a mismatch stops the run instead of producing a plausible wrong answer. + @test_throws BlockMismatchError create_parameters( + mesh, asm, (b1 = Poisson(f), b2 = Poisson(f)), one_props + ) + @test_throws BlockMismatchError create_parameters( + mesh, asm, one_physics, (b1 = one_props, b2 = one_props, nope = one_props) + ) + end +end