fix(fem): resolve an element surface to the face it names, not every node

src/ada/fem/shapes/shells.py

@@ -64,3 +64,47 @@
     ShellShapes.QUAD8: _QUAD_CORNER_FACES,
     ShellShapes.QUAD9: _QUAD_CORNER_FACES,
 }
+
+
+# ---------------------------------------------------------------------------
+# Abaqus shell edges, mid-side nodes included.
+#
+# ``shell_edges`` above is visualization topology (corner-to-corner only, and
+# the same list reused for the second-order shapes). A ``*Surface,
+# type=ELEMENT`` row naming ``E2`` needs exactly the nodes lying on that edge,
+# which for a second-order shell includes the edge's mid-side node.
+#
+# Ordering authority: Abaqus Analysis User's Guide, "Shell elements" ->
+# "Defining edge loads and surfaces on shells", where edge ``En`` runs between
+# corner nodes ``n`` and ``n+1`` (wrapping), read against the shell node
+# numbering in the same library section: the mid-side node of edge ``En`` is
+# node ``3+n`` on a 6-node triangle and ``4+n`` on an 8-node quadrilateral.
+# adapy's native ordering is Abaqus' here (see
+# ``node_order.NATIVE_MIDSIDE_EDGES``), so the 1-based numbers translate by
+# subtracting one.
+#
+# Indexed 0-based: entry ``i`` is Abaqus edge ``E(i+1)``.
+#
+# Note that ``SPOS`` / ``SNEG`` are *not* in these tables: they name a face of
+# the shell, i.e. the whole element, and every node of the element lies on it.
+# Callers handle those separately.
+
+_TRI_ABAQUS_EDGES = ((0, 1), (1, 2), (2, 0))
+_QUAD_ABAQUS_EDGES = ((0, 1), (1, 2), (2, 3), (3, 0))
+# 6-node triangle: node 4 on edge 1-2, 5 on 2-3, 6 on 3-1 (slots 3..5).
+_TRI6_ABAQUS_EDGES = ((0, 1, 3), (1, 2, 4), (2, 0, 5))
+# 8-node quadrilateral: node 5 on edge 1-2, 6 on 2-3, 7 on 3-4, 8 on 4-1
+# (slots 4..7). TRI7 / QUAD9 carry an extra *centre* node, which lies on no
+# edge, so they reuse the TRI6 / QUAD8 tables unchanged.
+_QUAD8_ABAQUS_EDGES = ((0, 1, 4), (1, 2, 5), (2, 3, 6), (3, 0, 7))
+
+#: ``{shape: (edge E1 slots, edge E2 slots, ...)}`` -- Abaqus shell edge
+#: numbering, mid-side nodes included.
+shell_abaqus_edges = {
+    ShellShapes.TRI: _TRI_ABAQUS_EDGES,
+    ShellShapes.TRI6: _TRI6_ABAQUS_EDGES,
+    ShellShapes.TRI7: _TRI6_ABAQUS_EDGES,
+    ShellShapes.QUAD: _QUAD_ABAQUS_EDGES,
+    ShellShapes.QUAD8: _QUAD8_ABAQUS_EDGES,
+    ShellShapes.QUAD9: _QUAD8_ABAQUS_EDGES,
+}

src/ada/fem/shapes/solids.py

@@ -132,3 +132,114 @@
     SolidShapes.WEDGE: _WEDGE_CORNER_FACES,
     SolidShapes.WEDGE15: _WEDGE_CORNER_FACES,
 }
+
+
+# ---------------------------------------------------------------------------
+# Abaqus element faces, mid-side nodes included.
+#
+# The ``solid_faces`` tables above are *visualization* topology: corner nodes
+# only, and for the wedge the quad faces are already split into triangles, so
+# the list is longer than the element's face count. Neither is usable for
+# resolving ``*Surface, type=ELEMENT`` -- a surface names one physical face
+# (``S3``) and needs exactly the nodes on it, mid-side nodes and all.
+#
+# Ordering authority: Abaqus Analysis User's Guide, "Three-dimensional solid
+# element library" -> the "Element faces" list published for each family
+# (C3D4/C3D10, C3D8/C3D20, C3D6/C3D15), read against the node numbering in the
+# same section. adapy's native node ordering *is* Abaqus' for these shapes
+# (see ``node_order.NATIVE_MIDSIDE_EDGES``), so the published 1-based numbers
+# translate to these slots by subtracting one. Within each face the nodes are
+# listed the way Abaqus lists them -- corners in face-winding order, then the
+# mid-side node of each of that face's edges in the same order -- so the tuples
+# double as a statement of the mid-side convention.
+#
+# Indexed 0-based: entry ``i`` is Abaqus face ``S(i+1)``.
+
+# C3D4 / C3D10 faces: S1 = 1-2-3, S2 = 1-4-2, S3 = 2-4-3, S4 = 3-4-1.
+_TETRA_ABAQUS_FACES = (
+    (0, 1, 2),
+    (0, 3, 1),
+    (1, 3, 2),
+    (2, 3, 0),
+)
+# C3D10 adds the mid-side node of each face edge: node 5 = edge 1-2, 6 = 2-3,
+# 7 = 1-3, 8 = 1-4, 9 = 2-4, 10 = 3-4 (slots 4..9).
+_TETRA10_ABAQUS_FACES = (
+    (0, 1, 2, 4, 5, 6),  # S1: edges 1-2, 2-3, 3-1
+    (0, 3, 1, 7, 8, 4),  # S2: edges 1-4, 4-2, 2-1
+    (1, 3, 2, 8, 9, 5),  # S3: edges 2-4, 4-3, 3-2
+    (2, 3, 0, 9, 7, 6),  # S4: edges 3-4, 4-1, 1-3
+)
+
+# C3D8 / C3D20 faces: S1 = 1-2-3-4, S2 = 5-8-7-6, S3 = 1-5-6-2,
+# S4 = 2-6-7-3, S5 = 3-7-8-4, S6 = 4-8-5-1.
+_HEX_ABAQUS_FACES = (
+    (0, 1, 2, 3),
+    (4, 7, 6, 5),
+    (0, 4, 5, 1),
+    (1, 5, 6, 2),
+    (2, 6, 7, 3),
+    (3, 7, 4, 0),
+)
+# C3D20 mid-side nodes: 9..12 on the 1-2-3-4 face edges, 13..16 on the
+# 5-6-7-8 face edges, 17..20 on the verticals 1-5, 2-6, 3-7, 4-8
+# (slots 8..19).
+_HEX20_ABAQUS_FACES = (
+    (0, 1, 2, 3, 8, 9, 10, 11),  # S1: edges 1-2, 2-3, 3-4, 4-1
+    (4, 7, 6, 5, 15, 14, 13, 12),  # S2: edges 5-8, 8-7, 7-6, 6-5
+    (0, 4, 5, 1, 16, 12, 17, 8),  # S3: edges 1-5, 5-6, 6-2, 2-1
+    (1, 5, 6, 2, 17, 13, 18, 9),  # S4: edges 2-6, 6-7, 7-3, 3-2
+    (2, 6, 7, 3, 18, 14, 19, 10),  # S5: edges 3-7, 7-8, 8-4, 4-3
+    (3, 7, 4, 0, 19, 15, 16, 11),  # S6: edges 4-8, 8-5, 5-1, 1-4
+)
+
+# C3D6 / C3D15 faces: S1 = 1-2-3, S2 = 4-5-6, S3 = 1-2-5-4,
+# S4 = 2-3-6-5, S5 = 3-1-4-6. Note only five faces -- the viz table above
+# has six entries because it triangulates the three quads.
+_WEDGE_ABAQUS_FACES = (
+    (0, 1, 2),
+    (3, 4, 5),
+    (0, 1, 4, 3),
+    (1, 2, 5, 4),
+    (2, 0, 3, 5),
+)
+# C3D15 mid-side nodes: 7..9 on the 1-2-3 face edges, 10..12 on the 4-5-6
+# face edges, 13..15 on the verticals 1-4, 2-5, 3-6 (slots 6..14).
+_WEDGE15_ABAQUS_FACES = (
+    (0, 1, 2, 6, 7, 8),  # S1: edges 1-2, 2-3, 3-1
+    (3, 4, 5, 9, 10, 11),  # S2: edges 4-5, 5-6, 6-4
+    (0, 1, 4, 3, 6, 13, 9, 12),  # S3: edges 1-2, 2-5, 5-4, 4-1
+    (1, 2, 5, 4, 7, 14, 10, 13),  # S4: edges 2-3, 3-6, 6-5, 5-2
+    (2, 0, 3, 5, 8, 12, 11, 14),  # S5: edges 3-1, 1-4, 4-6, 6-3
+)
+
+#: ``{shape: (face S1 slots, face S2 slots, ...)}`` -- Abaqus face numbering,
+#: mid-side nodes included. Deliberately incomplete: PYRAMID5 / PYRAMID13 and
+#: HEX27 are absent because Abaqus has no pyramid continuum element and no
+#: published C3D27 face-node numbering to copy, so there is no authority to
+#: follow for them. A caller reaching for a missing entry is meant to fail
+#: loudly rather than be handed a guess (or, worse, every node of the element).
+solid_abaqus_faces = {
+    SolidShapes.TETRA: _TETRA_ABAQUS_FACES,
+    SolidShapes.TETRA10: _TETRA10_ABAQUS_FACES,
+    SolidShapes.HEX8: _HEX_ABAQUS_FACES,
+    SolidShapes.HEX20: _HEX20_ABAQUS_FACES,
+    SolidShapes.WEDGE: _WEDGE_ABAQUS_FACES,
+    SolidShapes.WEDGE15: _WEDGE15_ABAQUS_FACES,
+}
+
+#: The same faces in the same order, corner slots only. A second-order element's faces
+#: are numbered exactly like its first-order counterpart's, so the corner slots of a
+#: C3D10 face are simply the C3D4 face of the same number -- no assumption needed about
+#: mid-side slots coming last inside a tuple. Used where a face has to be *identified*
+#: from a set of nodes (``surfaces.solid_abaqus_face_index``) rather than expanded into
+#: its nodes; matching on corners keeps that identification working for a caller that
+#: passes corner nodes only.
+solid_abaqus_face_corners = {
+    SolidShapes.TETRA: _TETRA_ABAQUS_FACES,
+    SolidShapes.TETRA10: _TETRA_ABAQUS_FACES,
+    SolidShapes.HEX8: _HEX_ABAQUS_FACES,
+    SolidShapes.HEX20: _HEX_ABAQUS_FACES,
+    SolidShapes.WEDGE: _WEDGE_ABAQUS_FACES,
+    SolidShapes.WEDGE15: _WEDGE_ABAQUS_FACES,
+}

src/ada/fem/surfaces.py

@@ -125,10 +125,21 @@ def create_surface_from_nodes(surface_name: str, nodes: List[Node], fem: "FEM",
 def get_surface_from_nodes_on_solid_elements(
     surface_name: str, all_el: List[Elem], nodes: List[Node], fem: "FEM", shell_positive: bool
 ) -> Surface:
+    """A solid ``Surface`` covering the faces of ``all_el`` that lie in ``nodes``.
+
+    The face number written into the set name and ``el_face_index`` is an **Abaqus** face
+    number, which is what every consumer of this surface expects -- including
+    :func:`side_node_indices`, which resolves it back into nodes. It used to come from
+    ``elem_has_parallel_face``, i.e. an index into the *visualisation* face sequence.
+    Those two agree for tetrahedra and hexahedra by coincidence, but not for a wedge: the
+    viz table triangulates the three quad faces, so it has six entries for a five-faced
+    element. A wedge could therefore come out labelled with the wrong face, or with an
+    "S6" that does not exist on that element at all.
+    """
     elements = []
     face_seq_indices = {}
     for el in all_el:
-        paralell_face_index = elem_has_parallel_face(el, nodes)
+        paralell_face_index = solid_abaqus_face_index(el, nodes)
         if paralell_face_index is None:
             continue
 
@@ -167,7 +178,36 @@ def get_surface_from_nodes_on_shell_elements(
     return Surface(surface_name, Surface.TYPES.ELEMENT, fs, el_face_index=side_name)
 
 
+def solid_abaqus_face_index(el: Elem, nodes: List[Node]) -> Union[int, None]:
+    """0-based Abaqus face number of the face of solid ``el`` whose corner nodes all lie
+    in ``nodes``, or ``None`` if no face does.
+
+    Matching is on corner nodes only, which is what the viz-table version did and what
+    keeps this working for a caller that collects corner nodes without the mid-side ones.
+    The number returned is the Abaqus one, so ``S{index + 1}`` names a face that really
+    exists on that element -- see :data:`ada.fem.shapes.solids.solid_abaqus_face_corners`.
+    """
+    from ada.fem.shapes.solids import solid_abaqus_face_corners
+
+    faces = solid_abaqus_face_corners.get(el.type)
+    if faces is None:
+        raise ValueError(
+            f"No Abaqus face numbering is defined for element type {el.type}. A surface on it "
+            f"cannot be given a face number; see solid_abaqus_faces for the types that are covered."
+        )
+    for i, slots in enumerate(faces):
+        if all(el.nodes[s] in nodes for s in slots):
+            return i
+    return None
+
+
 def elem_has_parallel_face(el: Elem, nodes: List[Node]):
+    """Whether any *visualisation* face of ``el`` lies wholly in ``nodes``, as its index.
+
+    Used by the shell path purely as a predicate -- the index is discarded there, since a
+    shell surface is named by side (``SPOS``/``SNEG``), not by face number. For solids use
+    :func:`solid_abaqus_face_index`, whose index is an Abaqus face number.
+    """
     for i, nid_refs in enumerate(el.shape.faces_seq):
         all_face_nodes_in_plane = True
         for nid in nid_refs:
@@ -180,36 +220,151 @@ def elem_has_parallel_face(el: Elem, nodes: List[Node]):
     return None
 
 
+#: Sides that cover the whole element rather than one of its faces. A shell's
+#: ``SPOS`` / ``SNEG`` are its two *faces*: every node of the element lies on
+#: each of them, so no filtering applies. The reader spells the same two as
+#: ``+1`` / ``-1`` when it has a shell elset in hand.
+WHOLE_ELEMENT_SIDES = frozenset({"SPOS", "SNEG"})
+
+
+def side_node_indices(el_type, side) -> Union[tuple, None]:
+    """Local node slots of ``el_type`` lying on surface side ``side``.
+
+    ``None`` means "the whole element" -- either no side was named at all, or the
+    side is a shell face (``SPOS`` / ``SNEG``), which every node of the element
+    lies on.
+
+    ``side`` is taken as the deck wrote it: ``"S3"`` / ``"E2"`` / ``"SPOS"``, or the
+    integer the Abaqus reader normalises a single-elset surface to (0-based face
+    index for a solid, ``+1`` / ``-1`` for a shell).
+
+    Raises on a side this has no map for, naming the element type and the side.
+    Falling back to every node of the element is exactly the bug this exists to
+    fix: a surface would silently cover the interior of its elements and, for a
+    tie or coupling, roughly twice the nodes the deck asked for.
+    """
+    from ada.fem.shapes.definitions import ShellShapes, SolidShapes
+    from ada.fem.shapes.shells import shell_abaqus_edges
+    from ada.fem.shapes.solids import solid_abaqus_faces
+
+    if side is None:
+        return None
+
+    if isinstance(side, str):
+        spec = side.strip().upper()
+        if spec == "" or spec in WHOLE_ELEMENT_SIDES:
+            return None
+        table, index = None, None
+        if len(spec) > 1 and spec[1:].isdigit():
+            index = int(spec[1:]) - 1
+            if spec[0] == "S":
+                table = solid_abaqus_faces.get(el_type)
+            elif spec[0] == "E":
+                table = shell_abaqus_edges.get(el_type)
+    else:
+        # The reader's normalised form. For a shell it is the +1/-1 SPOS/SNEG flag,
+        # not a face index; for a solid it is already 0-based.
+        if isinstance(el_type, ShellShapes):
+            return None
+        index = int(side)
+        table = solid_abaqus_faces.get(el_type) if isinstance(el_type, SolidShapes) else None
+
+    if table is None or index is None or not 0 <= index < len(table):
+        raise ValueError(
+            f'No Abaqus face/edge map for element type "{el_type}" side "{side}". '
+            "Add the element's face numbering to ada.fem.shapes (solid_abaqus_faces / "
+            "shell_abaqus_edges) rather than letting the surface fall back to every "
+            "node of every element."
+        )
+
+    return table[index]
+
+
+def _region_groups(region: Union["Surface", FemSet]):
+    """The region as ``(members, side)`` groups, one per set the region names.
+
+    Grouping rather than flattening keeps the side alongside the members it applies
+    to, and lets the caller resolve the side once per element *type* per group
+    instead of once per element.
+    """
+    if isinstance(region, FemSet):
+        # A set names no side -- it is the members themselves, whole.
+        return [(region.members, None)]
+
+    groups = []
+    fem_set = region.fem_set
+    fem_sets = fem_set if isinstance(fem_set, list) else [fem_set]
+    el_face_index = region.el_face_index
+    if isinstance(el_face_index, list):
+        if len(el_face_index) != len(fem_sets):
+            raise ValueError(
+                f'Surface "{region.name}" has {len(fem_sets)} FemSet(s) but {len(el_face_index)} el_face_index entries'
+            )
+        sides = el_face_index
+    else:
+        sides = [el_face_index] * len(fem_sets)
+
+    nodal = region.type == SurfTypes.NODE
+    for fs, side in zip(fem_sets, sides):
+        if fs is not None:
+            groups.append((fs.members, None if nodal else side))
+
+    id_refs = region.id_refs or []
+    if id_refs and region.parent is None:
+        raise ValueError(f'Surface "{region.name}" references sets by name but has no parent FEM')
+
+    # A multi-row surface (``_LIP_UNDERSIDE_S3, S3`` / ``_LIP_UNDERSIDE_S1, S1`` ...)
+    # is held entirely in ``id_refs``: the Abaqus reader leaves ``fem_set`` and
+    # ``el_face_index`` as None there and keeps each row's side label in the second
+    # slot of the row. For a NODE surface that slot is a weight factor instead.
+    for ref, ref_side in ((r[0], (None if nodal else r[1] if len(r) > 1 else None)) for r in id_refs):
+        if isinstance(ref, str):
+            sets = region.parent.sets
+            fs = sets.get_nset_from_name(ref) if nodal else sets.get_elset_from_name(ref)
+            groups.append((fs.members, ref_side))
+        else:
+            member = region.parent.nodes.from_id(ref) if nodal else region.parent.elements.from_id(ref)
+            groups.append(([member], ref_side))
+
+    return groups
+
+
 def surface_nodes(region: Union[Surface, FemSet]) -> List[Node]:
     """Unique nodes covered by a surface (or a plain set), in first-seen order.
 
     A ``Surface`` names its region through one or more ``FemSet`` objects, or -- when a
     deck listed several sets under a single surface -- through ``id_refs`` entries
     naming those sets or element / node ids outright. Callers that only need the nodes
     shouldn't have to know which of the three they got.
-    """
-    if isinstance(region, FemSet):
-        members = list(region.members)
-    else:
-        members = []
-        fem_set = region.fem_set
-        for fs in fem_set if isinstance(fem_set, list) else [fem_set]:
-            if fs is not None:
-                members += list(fs.members)
-
-        id_refs = region.id_refs or []
-        if id_refs and region.parent is None:
-            raise ValueError(f'Surface "{region.name}" references sets by name but has no parent FEM')
-        nodal = region.type == SurfTypes.NODE
-        for ref in (r[0] for r in id_refs):
-            if isinstance(ref, str):
-                sets = region.parent.sets
-                members += list((sets.get_nset_from_name(ref) if nodal else sets.get_elset_from_name(ref)).members)
-            else:
-                members.append(region.parent.nodes.from_id(ref) if nodal else region.parent.elements.from_id(ref))
 
+    For an element-based surface the side matters: ``_LIP_UNDERSIDE_S3, S3`` covers the
+    S3 face of each of those tetrahedra, not all ten of their nodes. Only the nodes on
+    the named face (or shell edge) come back, mid-side nodes included, per
+    :func:`side_node_indices`. A ``FemSet`` handed over directly names no side, and a
+    NODE-type surface's members are nodes already, so both keep their whole-membership
+    behaviour.
+    """
     nodes = {}
-    for m in members:
-        for n in getattr(m, "nodes", [m]):
-            nodes.setdefault(n.id, n)
+    for members, side in _region_groups(region):
+        # One small map per group, not per element: a set holds at most a couple of
+        # element types, and the surface side is constant across the set.
+        per_el_type = {}
+        for m in members:
+            m_nodes = getattr(m, "nodes", None)
+            if m_nodes is None:  # a Node, from a nodal set or an explicit node id
+                nodes.setdefault(m.id, m)
+                continue
+            if side is not None:
+                el_type = m.type
+                try:
+                    slots = per_el_type[el_type]
+                except KeyError:
+                    slots = per_el_type[el_type] = side_node_indices(el_type, side)
+                if slots is not None:
+                    for i in slots:
+                        n = m_nodes[i]
+                        nodes.setdefault(n.id, n)
+                    continue
+            for n in m_nodes:
+                nodes.setdefault(n.id, n)
     return list(nodes.values())