diff --git a/src/ada/fem/shapes/shells.py b/src/ada/fem/shapes/shells.py index 95b015ca6..36d6a8a61 100644 --- a/src/ada/fem/shapes/shells.py +++ b/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, +} diff --git a/src/ada/fem/shapes/solids.py b/src/ada/fem/shapes/solids.py index ff9f3860c..8f4a3a15e 100644 --- a/src/ada/fem/shapes/solids.py +++ b/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, +} diff --git a/src/ada/fem/surfaces.py b/src/ada/fem/surfaces.py index dc61f8b8c..49d0f2a6e 100644 --- a/src/ada/fem/surfaces.py +++ b/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,6 +220,115 @@ 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. @@ -187,29 +336,35 @@ def surface_nodes(region: Union[Surface, FemSet]) -> List[Node]: 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()) diff --git a/tests/core/fem/test_surface_nodes.py b/tests/core/fem/test_surface_nodes.py new file mode 100644 index 000000000..601f0d16c --- /dev/null +++ b/tests/core/fem/test_surface_nodes.py @@ -0,0 +1,278 @@ +"""``surface_nodes`` must honour the element face / edge a surface names. + +An element-based surface names a side -- ``S3`` on a tetrahedron, ``E1`` on a shell. +Returning every node of every member element instead covers the elements' interiors +and their adjoining faces, which for a tie or a coupling means roughly twice the +nodes the deck asked for, silently tied rigidly. + +The face and edge maps under test are Abaqus'; the node identities pinned below were +read off Abaqus' own numbering (see the comments in ``ada.fem.shapes.solids`` / +``shells``) and the C3D10 case matches what Abaqus itself selects on the project +deck this was written for. +""" + +from __future__ import annotations + +import pytest + +import ada +from ada.fem import FemSet, Surface +from ada.fem.surfaces import surface_nodes + + +def _ids(region) -> list[int]: + return [n.id for n in surface_nodes(region)] + + +def _tet10_fem() -> ada.FEM: + """One C3D10 whose node ids equal their 1-based Abaqus slot number. + + That makes the expected face node ids readable straight off Abaqus' + published face definitions: S1 = 1-2-3 + mid-side 5-6-7, and so on. + """ + # corners of a unit tet, then the six mid-side nodes in Abaqus order + corners = [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)] + mids = [(0.5, 0, 0), (0.5, 0.5, 0), (0, 0.5, 0), (0, 0, 0.5), (0.5, 0, 0.5), (0, 0.5, 0.5)] + nodes = [ada.Node(p, i + 1) for i, p in enumerate(corners + mids)] + fem = ada.FEM("tet", nodes=ada.api.containers.Nodes(nodes)) + fem.add_elem(ada.fem.Elem(1, nodes, "TETRA10")) + return fem + + +def _shell_fem() -> ada.FEM: + """One QUAD (node ids 1-4) and one TRI (node ids 3, 4, 5) sharing edge 3-4.""" + pts = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0), (0.5, 2, 0)] + nodes = [ada.Node(p, i + 1) for i, p in enumerate(pts)] + fem = ada.FEM("sh", nodes=ada.api.containers.Nodes(nodes)) + fem.add_elem(ada.fem.Elem(1, nodes[:4], "QUAD")) + fem.add_elem(ada.fem.Elem(2, [nodes[2], nodes[3], nodes[4]], "TRIANGLE")) + return fem + + +@pytest.mark.parametrize( + "side, expected", + [ + ("S1", [1, 2, 3, 5, 6, 7]), + ("S2", [1, 4, 2, 8, 9, 5]), + ("S3", [2, 4, 3, 9, 10, 6]), + ("S4", [3, 4, 1, 10, 8, 7]), + ], +) +def test_tet10_surface_returns_only_the_named_face(side, expected): + """Six of the ten nodes, and exactly the six Abaqus puts on that face.""" + fem = _tet10_fem() + fset = fem.add_set(FemSet(f"_s_{side}", [fem.elements.from_id(1)])) + surf = Surface(f"s_{side}", Surface.TYPES.ELEMENT, fset, el_face_index=side, parent=fem) + + assert _ids(surf) == expected + + +def test_tet10_reader_style_integer_face_index(): + """The Abaqus reader normalises ``S3`` to the 0-based index 2; same answer.""" + fem = _tet10_fem() + fset = fem.add_set(FemSet("_s", [fem.elements.from_id(1)])) + surf = Surface("s", Surface.TYPES.ELEMENT, fset, el_face_index=2, parent=fem) + + assert _ids(surf) == [2, 4, 3, 9, 10, 6] + + +def test_tet10_surface_is_not_every_node_of_the_element(): + """The defect this file exists for: 6 nodes on the face, not all 10.""" + fem = _tet10_fem() + fset = fem.add_set(FemSet("_s", [fem.elements.from_id(1)])) + surf = Surface("s", Surface.TYPES.ELEMENT, fset, el_face_index="S1", parent=fem) + + assert len(fem.elements.from_id(1).nodes) == 10 + assert len(surface_nodes(surf)) == 6 + + +def test_femset_passed_directly_keeps_every_node(): + """A plain set names no side, so it stays the whole membership.""" + fem = _tet10_fem() + fset = fem.add_set(FemSet("_all", [fem.elements.from_id(1)])) + + assert _ids(fset) == list(range(1, 11)) + + +def test_node_type_surface_is_unchanged(): + """A NODE surface's members are nodes already -- no face to honour.""" + nodes = [ada.Node((i, 0, 0), i + 1) for i in range(4)] + fem = ada.FEM("n", nodes=ada.api.containers.Nodes(nodes)) + fem.add_set(FemSet("a", nodes[:2], "nset")) + fem.add_set(FemSet("b", nodes[2:], "nset")) + surf = Surface("s", Surface.TYPES.NODE, None, id_refs=[("a", ""), ("b", "")], parent=fem) + + assert sorted(_ids(surf)) == [1, 2, 3, 4] + + +def test_shell_edge_surface_returns_only_the_edge(): + """``E2`` on a QUAD is corner nodes 2-3, not all four.""" + fem = _shell_fem() + fset = fem.add_set(FemSet("_e2", [fem.elements.from_id(1)])) + surf = Surface("s", Surface.TYPES.ELEMENT, fset, el_face_index="E2", parent=fem) + + assert _ids(surf) == [2, 3] + + +def test_shell_edge_surface_over_mixed_quad_and_tri(): + """One set, two element types: each looks its own edge table up.""" + fem = _shell_fem() + fset = fem.add_set(FemSet("_e2", [fem.elements.from_id(1), fem.elements.from_id(2)])) + surf = Surface("s", Surface.TYPES.ELEMENT, fset, el_face_index="E2", parent=fem) + + # QUAD E2 = nodes 2-3; TRI E2 = its 2nd and 3rd nodes, i.e. 4 and 5 + assert _ids(surf) == [2, 3, 4, 5] + + +def test_shell_face_sides_cover_the_whole_element(): + """SPOS / SNEG are faces of the shell -- every node lies on them.""" + fem = _shell_fem() + fset = fem.add_set(FemSet("_f", [fem.elements.from_id(1)])) + + for side in ("SPOS", "SNEG", 1, -1): + surf = Surface(f"s{side}", Surface.TYPES.ELEMENT, fset, el_face_index=side, parent=fem) + assert _ids(surf) == [1, 2, 3, 4] + + +def test_multi_face_surface_unions_its_faces_in_first_seen_order(): + """A surface holding a list of sets pairs each with its own face index.""" + fem = _tet10_fem() + el = fem.elements.from_id(1) + fs1 = fem.add_set(FemSet("_a", [el])) + fs2 = fem.add_set(FemSet("_b", [el])) + surf = Surface("s", Surface.TYPES.ELEMENT, [fs1, fs2], el_face_index=["S1", "S2"], parent=fem) + + # S1 = 1,2,3,5,6,7 then S2 = 1,4,2,8,9,5 contributing 4,8,9 as new + assert _ids(surf) == [1, 2, 3, 5, 6, 7, 4, 8, 9] + + +def test_unknown_face_index_raises_naming_type_and_side(): + """Never a silent fall back to every node -- that is the bug.""" + fem = _tet10_fem() + fset = fem.add_set(FemSet("_s", [fem.elements.from_id(1)])) + surf = Surface("s", Surface.TYPES.ELEMENT, fset, el_face_index="S9", parent=fem) + + with pytest.raises(ValueError, match="TETRA10.*S9"): + surface_nodes(surf) + + +def test_solid_face_side_on_a_shell_element_raises(): + """A shell has no ``S3``; refuse rather than return the whole element.""" + fem = _shell_fem() + fset = fem.add_set(FemSet("_s", [fem.elements.from_id(1)])) + surf = Surface("s", Surface.TYPES.ELEMENT, fset, el_face_index="S3", parent=fem) + + with pytest.raises(ValueError, match="QUAD.*S3"): + surface_nodes(surf) + + +ABAQUS_MULTI_FACE_DECK = """*Heading +** multi-row element surfaces, the shape the project deck writes +*Node +1, 0., 0., 0. +2, 1., 0., 0. +3, 0., 1., 0. +4, 0., 0., 1. +5, 0.5, 0., 0. +6, 0.5, 0.5, 0. +7, 0., 0.5, 0. +8, 0., 0., 0.5 +9, 0.5, 0., 0.5 +10, 0., 0.5, 0.5 +*Element, type=C3D10 +1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 +*Elset, elset=_TET_SURF_S1 +1, +*Elset, elset=_TET_SURF_S2 +1, +*Surface, type=ELEMENT, name=TET_SURF +_TET_SURF_S1, S1 +_TET_SURF_S2, S2 +*Solid Section, elset=_TET_SURF_S1, material=Steel +*Material, name=Steel +*Elastic +2.1e11, 0.3 +""" + + +def test_multi_row_surface_read_from_an_inp_honours_each_row(tmp_path): + """The path the project deck takes: the reader keeps the side labels in + ``id_refs`` (``fem_set`` and ``el_face_index`` stay None), so the face has to + be read off the row rather than off ``el_face_index``.""" + inp = tmp_path / "surf.inp" + inp.write_text(ABAQUS_MULTI_FACE_DECK) + + a = ada.from_fem(inp, "abaqus") + (part,) = [p for p in a.get_all_parts_in_assembly(True) if p.fem.surfaces] + surf = part.fem.surfaces["TET_SURF"] + + assert surf.fem_set is None and surf.el_face_index is None # it is all in id_refs + assert sorted(_ids(surf)) == [1, 2, 3, 4, 5, 6, 7, 8, 9] # S1 u S2; node 10 is not on either + + +def _wedge_fem() -> ada.FEM: + """One WEDGE, node ids 1..6 equal to their 1-based Abaqus slot. + + Triangles 1-2-3 at y=0 and 4-5-6 at y=1, so each Abaqus face is a distinct set of + nodes and the face a surface lands on is unambiguous. + """ + pts = [(0, 0, 0), (1, 0, 0), (0, 0, 1), (0, 1, 0), (1, 1, 0), (0, 1, 1)] + nodes = [ada.Node(p, i + 1) for i, p in enumerate(pts)] + fem = ada.FEM("wedge", nodes=ada.api.containers.Nodes(nodes)) + fem.add_elem(ada.fem.Elem(1, nodes, "WEDGE")) + return fem + + +def test_wedge_surface_is_named_with_an_abaqus_face_number(): + """The face number written for a solid surface must be Abaqus's, not the index of the + visualisation face sequence. + + The two agree for tets and hexes by coincidence. They cannot for a wedge: the viz + table triangulates the three quad faces, giving six entries for a five-faced element. + Taking the viz index labelled Abaqus face S2 (the 4-5-6 triangle) as "S2" only by + luck, and could name a face "S6", which no wedge has. + """ + from ada.fem.shapes.solids import solid_abaqus_faces + from ada.fem.surfaces import create_surface_from_nodes, solid_abaqus_face_index + + fem = _wedge_fem() + el = fem.elements.from_id(1) + + # every Abaqus face of this wedge is found, and found as itself + for index, slots in enumerate(solid_abaqus_faces[el.type]): + face_nodes = [el.nodes[s] for s in slots] + assert solid_abaqus_face_index(el, face_nodes) == index + + # ... and the quad face S4 (2-3-6-5) is named S4 in the surface it produces + s4_nodes = [el.nodes[s] for s in solid_abaqus_faces[el.type][3]] + surf = create_surface_from_nodes("quad_face", s4_nodes, fem) + assert surf.el_face_index == [3], "S4 is Abaqus face index 3" + assert any("_S4" in fs.name for fs in surf.fem_set) + # the face really does resolve back to those four nodes + assert {n.id for n in surface_nodes(surf)} == {n.id for n in s4_nodes} + + +def test_wedge_never_gets_a_face_number_it_does_not_have(): + """A wedge has five faces. The viz sequence has six entries, and its last one would + have been written as "S6".""" + from ada.fem.shapes.solids import solid_abaqus_faces + from ada.fem.surfaces import solid_abaqus_face_index + + fem = _wedge_fem() + el = fem.elements.from_id(1) + assert len(solid_abaqus_faces[el.type]) == 5 + assert len(el.shape.faces_seq) == 6, "the viz table still triangulates; that is why this matters" + for nodes in ([el.nodes[s] for s in slots] for slots in solid_abaqus_faces[el.type]): + assert solid_abaqus_face_index(el, nodes) < 5 + + +def test_tet10_surface_face_number_is_unchanged(): + """The tet path must be untouched by the wedge fix: its viz and Abaqus tables agree.""" + from ada.fem.shapes.solids import solid_abaqus_faces + from ada.fem.surfaces import solid_abaqus_face_index + + fem = _tet10_fem() + el = fem.elements.from_id(1) + for index, slots in enumerate(solid_abaqus_faces[el.type]): + corners = [el.nodes[s] for s in slots[:3]] + assert solid_abaqus_face_index(el, corners) == index