fix(fem): a Sesam settlement reaches Sestra, one BLDEP term per dof pair, and an *MPC with instance-qualified nodes

src/ada/fem/formats/abaqus/read/reader.py

@@ -1303,13 +1303,23 @@ def get_constraints_from_inp(bulk_str: str, fem: FEM) -> Dict[str, Constraint]:
                 mpc_dict[(f"{block_name}_{t.lower()}", t)] = mpc_dict.pop((block_name, t))
 
     def mpc_nodes(refs) -> list:
-        """The nodes an MPC names: a node id, or every node of a named set."""
+        """The nodes an MPC names: a node -- by bare id, or instance-qualified as ``p-1.4`` --
+        or every node of a named set.
+
+        The instance-qualified form is what Abaqus/CAE writes for an ``*MPC`` inside an
+        ``*Assembly``, and ``get_set_from_assembly`` resolves it to the *node*, not to a set:
+        ``BEAM, P-1.4, P-1.6`` reached ``ref.members`` and took the whole import down with
+        ``AttributeError: 'NodeProxy' object has no attribute 'members'``. The same deck written
+        flat, with bare ids, read fine, so only assembly decks were affected.
+        """
         nodes = []
         for ref in refs:
             if isinstance(ref, (int, np.integer)):
                 nodes.append(fem.nodes.from_id(int(ref)))
-            else:
+            elif isinstance(ref, FemSet):
                 nodes.extend(ref.members)
+            else:
+                nodes.append(ref)
         return nodes
 
     def get_mpc(mpc_name, mpc_type, mpc_values):

src/ada/fem/formats/sesam/write/not_held.py

@@ -73,11 +73,17 @@ def report_bcs(fems: Iterable[FEM], deck_fem: FEM | None = None) -> None:
                 rep.omitted(STAGE, "Bc", bc.name, "on a node outside the part's mesh, which is all the deck holds")
                 continue
             if any(m not in (None, 0, 0.0) for m in (bc.magnitudes or ())):
-                rep.approximated(
+                # Written in full now (BNBCD FIX code 2 plus a BNDISPL record), so this is a
+                # note, not an approximation. It stays because a settlement is loading in
+                # Sesam: it lands in one load case, which is a thing to know when reading the
+                # deck, and it is the one BC that does not read back as what was written --
+                # the Sesam reader has no BNDISPL card, so the value is lost on the way in.
+                rep.note(
                     STAGE,
                     "Bc",
                     bc.name,
-                    "a prescribed displacement is written as a fixed dof; the value is not written",
+                    "a prescribed displacement is written as BNBCD FIX code 2 plus a BNDISPL record "
+                    "in the first load case; Sesam carries a settlement as loading",
                     magnitudes=[m for m in bc.magnitudes],
                 )
             if bc.amplitude is not None:

src/ada/fem/formats/sesam/write/write_bcs.py

@@ -23,6 +23,11 @@
 or no metadata, falls through to the convention; if neither applies, no DOF gets code 4
 and nothing is logged. :func:`retained_dofs` is the entry point that applies that rule.
 
+A prescribed displacement (a settlement) takes two cards, not one: FIX code 2 on BNBCD says
+*which* DOFs are prescribed and a BNDISPL record in a load case says *by how much*. See
+:func:`prescribed_displacements` and :func:`bndispl_str`; neither card does anything without the
+other, measured against Sestra.
+
 Everything else in here is the BLDEP companion codes; see :func:`bnbcd_str`.
 """
 
@@ -45,10 +50,16 @@
 # every DOF carrying it becomes part of the superelement's external interface.
 FREE = 0
 FIXED = 1
-PRESCRIBED = 2  # never written: ada's Bc magnitudes are not carried into BNDISPL yet
+PRESCRIBED = 2  # a settlement; its value goes on a BNDISPL record, see bndispl_str
 DEPENDENT = 3
 RETAINED = 4
 
+#: BNDISPL's DTYPE. "Invalid DTYPE on BNDISPL card. Supported values are 1 (displacement) and
+#: 3 (acceleration)" -- Sestra V11.3-00's own diagnostic, in ``Bin/DataAccess.dll``. A ``Bc``
+#: magnitude is a displacement; DTYPE 3 written instead leaves a static run with no
+#: displacement result at all (measured, see :func:`bndispl_str`).
+DTYPE_DISPLACEMENT = 1
+
 #: ``to_fem(..., metadata={RETAINED_KEY: {"<nset name>": [1, 2, 3, 4, 5, 6]}})``
 RETAINED_KEY = "sesam_retained_dofs"
 
@@ -247,19 +258,99 @@ def _fixed_dofs(bc) -> set[int]:
     return {int(d) for d in dofs if d is not None and not isinstance(d, str)}
 
 
+def prescribed_displacements(fems: Sequence[FEM]) -> dict[int, dict[int, float]]:
+    """``{node id: {dof: value}}`` -- the settlements that become BNDISPL records.
+
+    A ``Bc`` carrying ``magnitudes`` is a *prescribed* displacement, not a support: the Abaqus
+    data line ``push, 3, 3, -0.02`` moves that node 20 mm down. It used to be written as an
+    ordinary clamp (FIX code 1) with the value dropped, so a settlement silently became a rigid
+    support -- a different structure, and one whose reaction forces are the ones the engineer
+    is usually after.
+
+    Only a nonzero magnitude counts. A prescribed zero *is* a fixed support, so writing it as
+    code 2 plus a BNDISPL of 0.0 would change the text of every deck that has one without
+    changing the model it describes.
+
+    ``fems`` is the same sequence :func:`bnbcd_str` takes -- the writer passes the BCs BNBCD can
+    hold (``not_held.bc_is_held``), so a velocity BC's magnitude cannot arrive here as a
+    displacement. Two BCs prescribing one DOF differently leave the later one's value, which is
+    how the fixed DOFs of two BCs already merge.
+    """
+    out: dict[int, dict[int, float]] = {}
+    for fem in fems:
+        if fem is None:
+            continue
+        for bc in fem.bcs:
+            if bc.fem_set.type != "nset":
+                continue
+            for dof, magnitude in zip(bc.dofs, bc.magnitudes or ()):
+                if isinstance(dof, str) or magnitude is None or float(magnitude) == 0.0:
+                    continue
+                for mem in bc.fem_set.members:
+                    out.setdefault(int(mem.id), {})[int(dof)] = float(magnitude)
+    return out
+
+
+def bndispl_str(prescribed: dict[int, dict[int, float]], ndofs: NodeDofs | None = None, lid: int = 1) -> str:
+    """The BNDISPL block: one record per node with a prescribed displacement, node id order.
+
+    ``BNDISPL  LLC  DTYPE  COMPLX  0`` then ``NODENO  NDOF`` and NDOF values -- the same
+    four-field FORTRAN rows as BNLOAD (``write_loads.load_force``). ``COMPLX`` 0: no phase
+    shift, so no imaginary half follows ("BNDISPL card is too short for complex loads and dof
+    count given in NDOF field", Sestra V11.3-00's own diagnostic).
+
+    The record lives in a **load case** (``LLC``), which is the part of this that is not obvious:
+    a prescribed displacement is loading in Sesam, not boundary data. Measured on a 1 m IPE300
+    cantilever, tip dz = -0.01 prescribed, Sestra V11.3-00 solving it:
+
+    * FIX code 2 **and** this record: the SIN gives the tip Z = -0.01 exactly, and its rotation
+      free (RY = 6.907e-3), which is the settlement asked for.
+    * FIX code 2, no BNDISPL: "WARNING ... No load is specified", no displacement result at all.
+    * this record, FIX code 1 instead of 2: the value is ignored, tip Z = 0.0.
+
+    So the two cards are one statement and neither half is optional. ``lid`` is the load case
+    number the caller is writing (``write_loads.step_loads_str`` owns that numbering).
+    """
+    from .writer import ALL_SIX_DOF
+
+    if ndofs is None:
+        ndofs = ALL_SIX_DOF
+
+    out = ""
+    for nid in sorted(prescribed or {}):
+        ndof = ndofs.ndof(nid)
+        beyond = sorted(dof for dof in prescribed[nid] if dof > ndof)
+        if beyond:
+            raise ValueError(
+                f"sesam writer: a displacement is prescribed on dof(s) {beyond} of node {nid}, which has "
+                f"{ndof} dofs (NDOF={ndof}). A node attached only to solid elements has no rotational "
+                "dofs for a prescribed rotation to act on."
+            )
+        vals = [prescribed[nid].get(dof, 0.0) for dof in range(1, ndof + 1)]
+        out += write_ff(
+            "BNDISPL",
+            [(lid, DTYPE_DISPLACEMENT, 0, 0), (nid, ndof) + tuple(vals[:2]), tuple(vals[2:ndof])],
+        )
+    return out
+
+
 def bnbcd_str(
     fems: Sequence[FEM],
     lin_deps: Sequence[BldepRecord] = (),
     retained: dict[int, tuple[int, ...]] | None = None,
     ndofs: NodeDofs | None = None,
+    prescribed: dict[int, dict[int, float]] | None = None,
 ) -> str:
     """The BNBCD block: one record per node that has anything to say about its DOFs.
 
     Three sources are merged into the six FIX codes of each node:
 
-    * ``fem.bcs`` -> code 1 on every DOF the boundary condition names. Two BCs touching
-      the same node merge into one record holding the union of their DOFs (the writer
-      used to emit one record per BC, which left the same node declared twice).
+    * ``fem.bcs`` -> code 1 on every DOF the boundary condition names, or code 2 where
+      ``prescribed`` (``{node id: {dof: value}}`` from :func:`prescribed_displacements`)
+      carries a value for it -- a settlement, whose value goes on the BNDISPL record
+      :func:`bndispl_str` writes. Two BCs touching the same node merge into one record
+      holding the union of their DOFs (the writer used to emit one record per BC, which
+      left the same node declared twice).
     * ``lin_deps`` (the BLDEP records, see ``write_constraints.BldepRecord``) -> code 3
       on each dependent DOF of the dependent node, and an explicit all-free record for
       each independent (master) node. BLDEP is only valid alongside these codes, manual
@@ -322,12 +413,13 @@ def check_dof(nid: int, dof: int, what: str) -> None:
             fixed_dofs = _fixed_dofs(bc)
             for mem in bc.fem_set.members:
                 node = node_codes(mem.id)
+                settled = (prescribed or {}).get(int(mem.id), {})
                 for dof in range(1, 7):
-                    # As before: a named DOF is written as fixed regardless of its
-                    # magnitude. Prescribed values (code 2) are not emitted by ada.
+                    # A named DOF is fixed unless a magnitude was given for it, in which case
+                    # it is prescribed and BNDISPL carries the value.
                     if dof in fixed_dofs:
                         check_dof(mem.id, dof, f'boundary condition "{bc.name}"')
-                        node[dof - 1] = FIXED
+                        node[dof - 1] = PRESCRIBED if dof in settled else FIXED
 
     chained_master_warned = False
     for record in lin_deps:

src/ada/fem/formats/sesam/write/write_constraints.py

@@ -65,7 +65,9 @@ def bldep_records(fem: FEM, ndofs: NodeDofs | None = None) -> list[BldepRecord]:
     if ndofs is None:
         ndofs = ALL_SIX_DOF
     rep = report()
-    records: list[BldepRecord] = []
+    # Each record is carried with the name of the constraint that produced it: ``_merged`` has to
+    # say which two constraints disagree when they claim one (slave_dof, master_dof).
+    owned: list[tuple[str, BldepRecord]] = []
     claimed: dict[tuple[int, int], str] = {}
     for constraint in fem.constraints.values():
         # A rigid body links every slave node rigidly to the master (reference) node — the
@@ -90,34 +92,93 @@ def bldep_records(fem: FEM, ndofs: NodeDofs | None = None) -> list[BldepRecord]:
         else:
             rep.omitted(STAGE, "Constraint", constraint.name, f'a "{constraint.type}" constraint has no Sesam form')
             continue
-        records += new
-        for key in {(r.slave, dof) for r in new for dof in r.slave_dofs}:
+        owned += [(constraint.name, r) for r in new]
+        for key in sorted({(r.slave, dof) for r in new for dof in r.slave_dofs}):
             first = claimed.setdefault(key, constraint.name)
             if first != constraint.name:
                 rep.suspect(
                     STAGE,
                     "Constraint",
                     constraint.name,
-                    "makes a dof dependent that another constraint already does; Sesam sums linear dependencies",
+                    "makes a dof dependent that another constraint already does; nothing in the "
+                    "deck says which of the two relations is meant",
                     node=key[0],
                     dof=key[1],
                     other=first,
                 )
 
-    return _merged(records)
+    return _merged(owned)
 
 
-def _merged(records: list[BldepRecord]) -> list[BldepRecord]:
-    """One record per (slave, master) pair, in first-seen order.
+#: How close two betas have to be to count as the same dependency declared twice rather than two
+#: different ones. Duplicates come from the same arithmetic (the same ``*Equation`` read twice,
+#: the same lever arm through ``LinDep``) and agree bit for bit; the window is here only so that
+#: float noise is not read as a deliberately different coefficient.
+_BETA_SAME = 1e-12
+
+
+def _merged(owned: list[tuple[str, BldepRecord]]) -> list[BldepRecord]:
+    """One record per (slave, master) pair, and one term per (slave_dof, master_dof), in
+    first-seen order. ``owned`` pairs each record with the name of the constraint it came from.
 
     The manual (BLDEP, section 7.2.14): "The same combination of SLAVE and MASTER may occur
     only once." Two equations relating the same pair of nodes -- or an equation with two
     terms on one independent node -- therefore share one record holding all their terms.
+
+    A ``(slave_dof, master_dof)`` may not repeat *inside* that record either, and used to: the
+    terms were concatenated, and Sestra adds the betas of a repeated pair together. A deck
+    declaring ``u(1,1) = 1.0 * u(2,1)`` twice came out as ``NDDOF 1 NDEP 2`` with beta 1.0
+    twice, which is ``u(1,1) = 2.0 * u(2,1)`` -- the dependency doubled, with nothing in the
+    written deck to show it had been meant once. What produced the repeat decides what to do
+    with it:
+
+    * the same beta twice is one dependency declared twice (the same ``*Equation`` pasted in
+      again, two Abaqus keywords the reader turned into the same relation). Written once, and
+      noted -- nothing is lost.
+    * two different betas are two incompatible statements about one dof: a ``*Tie`` and an
+      ``*Equation`` over the same node pair, which is a modelling error in the source deck.
+      Summing them is not what either constraint says and neither is averaging, so the later
+      one is refused by name in the report and the first-declared relation is what is written.
     """
-    merged: dict[tuple[int, int], list] = {}
-    for r in records:
-        merged.setdefault((r.slave, r.master), []).extend(r.terms)
-    return [BldepRecord(s, m, tuple(terms)) for (s, m), terms in merged.items()]
+    rep = report()
+    merged: dict[tuple[int, int], dict[tuple[int, int], float]] = {}
+    owners: dict[tuple[int, int, int, int], str] = {}
+    for name, r in owned:
+        terms = merged.setdefault((r.slave, r.master), {})
+        for s_dof, m_dof, beta in r.terms:
+            first = owners.setdefault((r.slave, r.master, s_dof, m_dof), name)
+            if (s_dof, m_dof) not in terms:
+                terms[(s_dof, m_dof)] = beta
+                continue
+            kept = terms[(s_dof, m_dof)]
+            if abs(beta - kept) <= _BETA_SAME * max(1.0, abs(beta), abs(kept)):
+                rep.note(
+                    STAGE,
+                    "Constraint",
+                    name,
+                    "declares a linear dependency another constraint already declares; it is "
+                    "written once, as Sestra would otherwise add the two together",
+                    node=r.slave,
+                    master=r.master,
+                    dof=s_dof,
+                    beta=beta,
+                    other=first,
+                )
+            else:
+                rep.omitted(
+                    STAGE,
+                    "Constraint",
+                    name,
+                    "gives a dof a different dependency on the same master dof than another "
+                    "constraint already does; BLDEP holds one, so the first one is kept",
+                    node=r.slave,
+                    master=r.master,
+                    dof=s_dof,
+                    beta=beta,
+                    kept=kept,
+                    other=first,
+                )
+    return [BldepRecord(s, m, tuple((sd, md, b) for (sd, md), b in terms.items())) for (s, m), terms in merged.items()]
 
 
 def _report_coupling(constraint: Constraint, records: list[BldepRecord], ndofs: NodeDofs) -> None:

src/ada/fem/formats/sesam/write/write_loads.py

@@ -15,37 +15,75 @@
     from .writer import NodeDofs
 
 
-def loads_str(fem: FEM, ndofs: NodeDofs | None = None) -> str:
+#: The load case a step's loads go into when the step declares none of its own.
+DEFAULT_CASE = "LC1"
+
+
+def loads_str(fem: FEM, ndofs: NodeDofs | None = None, prescribed: dict[int, dict[int, float]] | None = None) -> str:
     """The load block of ``fem``'s first step.
 
     ``ndofs`` (:class:`writer.NodeDofs`) is threaded down to BNLOAD, which declares an
     NDOF of its own and must agree with GNODE. Left at ``None`` it is derived from ``fem``.
     """
     from .writer import node_dofs
 
-    if len(fem.steps) == 0:
-        return ""
-    return step_loads_str(fem.steps[0], node_dofs(fem) if ndofs is None else ndofs)
+    step = fem.steps[0] if len(fem.steps) > 0 else None
+    return step_loads_str(step, node_dofs(fem) if ndofs is None else ndofs, prescribed)
 
 
-def step_loads_str(step: Step | None, ndofs: NodeDofs | None = None) -> str:
+def _case_str(lid: int, name: str) -> str:
+    return write_ff("TDLOAD", [(4, lid, 100 + len(name), 0), (name,)])
+
+
+def step_loads_str(
+    step: Step | None,
+    ndofs: NodeDofs | None = None,
+    prescribed: dict[int, dict[int, float]] | None = None,
+) -> str:
     """The load block of one step: its load cases, or all its loads as one case ``LC1``.
 
     ``to_fem`` passes the part FEM's dof counts even when the step lives on the assembly,
     because that is where the nodes a load names actually live.
+
+    ``prescribed`` (``{node id: {dof: value}}``, from ``write_bcs.prescribed_displacements``)
+    are the settlements, and they are written *here* rather than with the boundary conditions
+    because in Sesam a prescribed displacement is loading: its BNDISPL record declares an LLC.
+    A model whose only loading is a settlement therefore still needs a load case -- with FIX
+    code 2 and no load case at all Sestra V11.3-00 warns "No load is specified" and writes no
+    displacement result -- so one is opened for it here.
     """
+    from .write_bcs import bndispl_str
+
+    prescribed = prescribed or {}
     if step is None or len(step.loads) == 0:
-        return ""
+        if not prescribed:
+            return ""
+        return _case_str(1, DEFAULT_CASE) + bndispl_str(prescribed, ndofs, 1)
 
     if len(step.load_cases.keys()) > 0:
         cases = [(lc.name, lc.loads or []) for lc in step.load_cases.values()]
     else:
-        cases = [("LC1", step.loads)]
+        cases = [(DEFAULT_CASE, step.loads)]
 
     out_str = ""
     for lid, (lc_name, loads) in enumerate(cases, start=1):
-        out_str += write_ff("TDLOAD", [(4, lid, 100 + len(lc_name), 0), (lc_name,)])
+        out_str += _case_str(lid, lc_name)
         out_str += case_loads_str(loads, lid, ndofs)
+        if lid == 1:
+            out_str += bndispl_str(prescribed, ndofs, lid)
+    if prescribed and len(cases) > 1:
+        # A Bc belongs to no load case, so which case a settlement acts in is not something the
+        # model says. Sestra solves each case on its own, so putting it in all of them would make
+        # every case a settlement case; the first one is the choice, said out loud.
+        report().approximated(
+            STAGE,
+            "BNDISPL",
+            cases[0][0],
+            "the prescribed displacements are written into the first load case only; each Sesam "
+            "load case has its own BNDISPL records and a Bc belongs to no load case",
+            n_nodes=len(prescribed),
+            n_cases=len(cases),
+        )
     return out_str