BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
Thin-shell virtual domains and edge allocation

Date: 2026-06-16 Purpose: Explain what a virtual domain is in BELFEM's thin-shell formulation, why its volume elements exist only in the output mesh, and the rule that decides whether such elements carry edges — which governs how the periodic rebuild and other edge-consuming passes must treat them. Module: src/fem/maxwell


What a virtual domain is

A virtual domain is a set of element blocks whose volumes do not exist in the original input geometry. They are generated by BELFEM, not read from the mesh.

For a thin shell (a superconducting tape, film, or other thin layer), the input geometry contains only a surface. BELFEM extrudes that surface into one or more volume layers using the thicknesses provided in input.conf. The resulting extruded volume blocks are the virtual domain of the thin shell.

These virtual blocks are added to the output .exo mesh for visualization purposes — so the extruded shell can be seen and post-processed like a real volume. Because they are extruded from a surface that sits inside the model, the virtual blocks geometrically intersect the surrounding air elements. That overlap is expected and correct — it is exactly how the thin-shell formulation is meant to work (the shell is a degenerate-thickness object represented on a surface; the extrusion is a visualization/representation device, not a true cut-out of the air).


Edge allocation rule (critical)

Whether a thin-shell (virtual) block's elements receive edges depends on the block's domain type:

thin-shell block domain type gets edges? why
Conductor yes carries the H-field edge (h) DOFs of the H-φ formulation
Ferro no non-conducting; no edge DOFs
Void no non-conducting; no edge DOFs
Buffer (non-conducting) no non-conducting; no edge DOFs

So a virtual block's elements have edge containers (Element::has_edges() true) only when the block is a conductor. Ferro, void, and non-conducting buffer virtual elements are created without edges.

This is consistent with the H-φ formulation: only conducting regions carry the edge-based H field; non-conducting regions are described by the scalar potential φ on nodes and need no edges.


Consequence for edge-consuming passes (e.g. the periodic rebuild)

Any pass that iterates a facet's edges via Facet::edge() — which delegates to the wrapper element and is guarded by a debug mHaveEdges assert (cl_ElementTemplate.hpp, edge()) — must not touch facets whose wrapper element belongs to a non-conducting virtual block.

The concrete case that surfaces this: PeriodicityFactory::collect_edges() iterates the edges of the selected periodic seam facets. When a thin-shell virtual block lies on a periodic plane, its seam facets are selected geometrically, but Facet::edge() then asserts (Edges for element N on block B have not been allocated).

Important distinction (wrapper facet vs volume element). The edge-allocation rule above governs the volume elements of a virtual block. The generated periodic wrapper facets are a separate matter: ThinShellFactory::create_periodic_sideset() builds them with only set_master()/set_sideset_id() and never links them to any edge container — so a thin-shell periodic wrapper has no edges regardless of whether its backing volume is a conductor. (Facet::edge() delegates to the wrapper element, not the master volume element.) This is compounded by ordering: the full edge build and first periodicity run in CutFactory before the thin shell exists; create_edges_and_faces_on_mesh() later rebuilds conductor volume edges only (no sideset wrappers); then the Maxwell periodic update re-collects edges from the now-present thin-shell wrappers.

Therefore: a bare ! element()->has_edges() skip in collect_edges() would be only a crash guard — it would also skip conductor thin-shell periodic wrappers and silently drop their periodic edge constraints.

How this is actually solved: the master-facet fallback. PeriodicityFactory::collect_edges() (src/mesh/cl_Mesh_PeriodicityFactory.cpp:600) does not require the wrapper to own edges. It tries the wrapper element first, and when that has no edge container it reconstructs the edge set from the master volume element via master()->get_edges_of_facet( index_on_master() ) (:624-635). So a conductor wrapper's periodic edges are contributed even though create_periodic_sideset() never allocates its edge container.

Two things make the reconstruction safe, and both are load-bearing:

  • Enumeration and intrinsic direction may differ from the partner side. That is tolerated rather than prevented: match_edges() (:886) pairs edges by their endpoint originals and aligns direction by node swap, so neither the ordering nor the sign of the master-derived edges has to agree across the seam.
  • The master is always a TET4. Because the shell normal is never coplanar with a periodic plane, every 3D periodic seam facet is a TRI3 backed by a TET4 (see src/mesh/doc/thin_shell_geometry_and_periodicity.md §5.1). This matters because TET4::get_edges_of_facet is defined for all four faces, whereas PENTA6TS::get_edges_of_facet covers only the triangular caps — the QUAD4 laterals deliberately carry no through-thickness edges in the reduced TS edge model. Were a PENTA lateral ever to become a seam facet, the fallback would have no edges to source.

An alternative fix — linking each generated wrapper facet to the thin-shell layer edges at creation time — was considered and not taken. It is a thin-shell-aware operation (the side facets of PENTA*TS elements have their own edge enumeration), and the master-facet fallback achieves the same result inside the one function that consumes the edges.