Give element blocks one canonical order, and match physics by name

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 <cmhamel32@gmail.com> and contributors"]
 
 [deps]

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

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"]

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

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

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)