diff --git a/web/package.json b/web/package.json index a77798d..76ea24e 100644 --- a/web/package.json +++ b/web/package.json @@ -17,7 +17,7 @@ "examples:generate": "bun ../examples/web-examples/generate.mjs", "examples:generate-crane-shell": "bun ../examples/web-examples/generate-crane-shell.mjs", "examples:generate-plate-hole-shell": "bun ../examples/web-examples/generate-plate-hole-shell.mjs", - "test": "bun ../examples/validation/run.mjs && bun ../examples/validation/dof-constraint.test.mjs && bun ../examples/validation/prescribed-displacement.test.mjs && bun tests/test_wall_bracket.mjs 20.0 && bun tests/test_multibody.mjs && bun tests/test_tie.mjs && bun tests/test_coupling.mjs && bun tests/test_face_pick.mjs && bun tests/test_edge_pick.mjs && bun tests/test_moment_load.mjs && bun tests/test_shellize_mpc.mjs && bun tests/test_body_highlight.mjs && bun tests/test_midsurface_junction.mjs && bun tests/test_coupled_materials.mjs && playwright test", + "test": "bun ../examples/validation/run.mjs && bun ../examples/validation/dof-constraint.test.mjs && bun ../examples/validation/prescribed-displacement.test.mjs && bun tests/test_wall_bracket.mjs 20.0 && bun tests/test_multibody.mjs && bun tests/test_tie.mjs && bun tests/test_coupling.mjs && bun tests/test_reference_point.mjs && bun tests/test_face_pick.mjs && bun tests/test_edge_pick.mjs && bun tests/test_moment_load.mjs && bun tests/test_shellize_mpc.mjs && bun tests/test_body_highlight.mjs && bun tests/test_midsurface_junction.mjs && bun tests/test_coupled_materials.mjs && playwright test", "test:coverage": "rm -rf .nyc_output coverage && COVERAGE=1 playwright test && bun run coverage:report", "coverage:report": "nyc report && bun scripts/coverage-dead-code.ts", "test:ui": "playwright test --ui", diff --git a/web/src/components/panel/BcLoadFormControls.tsx b/web/src/components/panel/BcLoadFormControls.tsx index 37aa66e..556b4ff 100644 --- a/web/src/components/panel/BcLoadFormControls.tsx +++ b/web/src/components/panel/BcLoadFormControls.tsx @@ -2,8 +2,11 @@ // SPDX-License-Identifier: AGPL-3.0-or-later import type { + CouplingKind, FaceSelection, LoadKind, + PickGeometry, + ReferencePointOption, TieExtent, } from "../../store/modelStore"; import { @@ -14,33 +17,51 @@ import { } from "./bcFormUtils"; import styles from "./LeftPanel.module.css"; -// Face / Edge picker for shell models: a flat shell sheet is one face-pick -// region, so grabbing its rim (an edge load or supported-edge BC) needs edge -// picking (#386). Offered only when the model has CTRIA3 elements. +const GEOMETRY_LABEL: Record = { + face: "Face", + edge: "Edge", + point: "Point", +}; + +// What a click selects. `options` is explicit because not every geometry suits +// every model or every condition: edge picking only means something on a shell +// (a closed solid boundary has no boundary polyline to walk), and a point pick +// is not offered where a single node would be the wrong thing to define. +// Rendered only when there is a genuine choice to make. export function PickGeometryToggle({ value, + options, onChange, }: { - value: "face" | "edge"; - onChange(geometry: "face" | "edge"): void; + value: PickGeometry; + options: PickGeometry[]; + onChange(geometry: PickGeometry): void; }) { + if (options.length < 2) return null; return (
- {(["face", "edge"] as const).map((geometry) => ( + {options.map((geometry) => ( ))}
); } +const PICK_HINT: Record = { + face: "Click a face in the 3D viewport", + edge: "Click near a mesh edge in the 3D viewport", + point: "Click near a mesh node in the 3D viewport", +}; + export function PickedFaceList({ faces, onRemove, @@ -48,16 +69,10 @@ export function PickedFaceList({ }: { faces: FaceSelection[]; onRemove(index: number): void; - geometry?: "face" | "edge"; + geometry?: PickGeometry; }) { if (faces.length === 0) { - return ( -
- {geometry === "edge" - ? "Click near a mesh edge in the 3D viewport" - : "Click a face in the 3D viewport"} -
- ); + return
{PICK_HINT[geometry]}
; } return (
@@ -229,3 +244,92 @@ export function DofCheckboxes({
); } + +// How a surface-to-point coupling ties its surface to its reference point. The +// two kinds are not interchangeable — see lib/coupling.ts — so the difference +// is spelled out under the select rather than left to the two words. +export function CouplingKindSelect({ + value, + onChange, +}: { + value: CouplingKind; + onChange(kind: CouplingKind): void; +}) { + return ( + <> +
+ Kind + +
+
+ {value === "distributing" + ? "the surface stays flexible and the point follows it — the point can be loaded, but not fixed" + : "the surface follows the point rigidly — the point can be fixed, loaded or coupled onward"} +
+ + ); +} + +// Where the reference point sits: one of the positions derived from the picked +// surface (its centre, and the two ends of its axis when it is a cylinder), or +// coordinates typed in. Choosing a derived position fills the coordinate boxes, +// which stay editable — the derived positions are a starting point, not a +// constraint (KOF-208). +export function ReferencePointInputs({ + options, + coords, + onCoordChange, + onPickOption, +}: { + options: ReferencePointOption[]; + coords: [string, string, string]; + onCoordChange(index: number, value: string): void; + onPickOption(option: ReferencePointOption): void; +}) { + return ( + <> + {options.length > 0 && ( +
+ Place at + +
+ )} + {["X", "Y", "Z"].map((axis, i) => ( +
+ {axis} (mm) + onCoordChange(i, e.target.value)} + /> +
+ ))} + + ); +} diff --git a/web/src/components/panel/BcSection.tsx b/web/src/components/panel/BcSection.tsx index ad37035..3c37bdb 100644 --- a/web/src/components/panel/BcSection.tsx +++ b/web/src/components/panel/BcSection.tsx @@ -26,6 +26,7 @@ export function BcSection({ onError }: { onError(msg: string | null): void }) { const hasShells = useModelStore((s) => s.elements.some((el) => el.type === "CTRIA3"), ); + const couplingGroups = useModelStore((s) => s.couplingGroups); const createBcGroup = useModelStore((s) => s.createBcGroup); const addFaceToBcGroup = useModelStore((s) => s.addFaceToBcGroup); const removeFaceFromBcGroup = useModelStore((s) => s.removeFaceFromBcGroup); @@ -51,6 +52,15 @@ export function BcSection({ onError }: { onError(msg: string | null): void }) { false, ]); const [bcValue, setBcValue] = useState("0"); + // Whether the pick has landed on a coupling's reference point — clicked in + // the viewport like anything else. Such a point carries six real DOFs + // (shell_core gives a coupling reference its rotations), so the rotational + // boxes become meaningful: fixing them is how a bolted or welded connection + // is restrained without clamping every node of the bore. + const refPoints = new Set(couplingGroups.map((c) => c.refNodeId)); + const pickedReferencePoint = allPickedFaces.some((face) => + face.nodeIds.every((id) => refPoints.has(id)), + ); const targetBcGroup = pickTargetGroupId !== null @@ -63,7 +73,7 @@ export function BcSection({ onError }: { onError(msg: string | null): void }) { allPickedFaces, // eslint-disable-next-line kofem/no-silent-fallback -- numbering offset for the new entries; a pick with no target group starts a fresh group, which has 0 faces targetBcGroup?.faces.length ?? 0, - pickGeometry === "edge" ? "Edge" : "Face", + pickGeometry, ); if (targetBcGroup) { for (const faceEntry of faceEntries) { @@ -109,12 +119,16 @@ export function BcSection({ onError }: { onError(msg: string | null): void }) { - {hasShells && ( - - )} + {/* Point picking is offered so a coupling's REFERENCE POINT can be + fixed — clicked in the viewport like any other selection, which is + the only way to restrain a bolted hole without clamping every node + of its bore. It also allows a single mesh node, which is a + legitimate way to remove a rigid-body mode. */} + 0 && !targetBcGroup && ( <> + {/* A reference point carries six DOFs whichever elements the + model has, so Rx/Ry/Rz appear as soon as one is picked. */} setCheckedDofs((prev) => prev.map((checked, i) => diff --git a/web/src/components/panel/BoundaryConditionsPanel.tsx b/web/src/components/panel/BoundaryConditionsPanel.tsx index c52a943..23ba443 100644 --- a/web/src/components/panel/BoundaryConditionsPanel.tsx +++ b/web/src/components/panel/BoundaryConditionsPanel.tsx @@ -6,6 +6,7 @@ import { useModelStore } from "../../store/modelStore"; import { BcSection } from "./BcSection"; import { LoadSection } from "./LoadSection"; import { TieSection } from "./TieSection"; +import { CouplingSection } from "./CouplingSection"; import styles from "./LeftPanel.module.css"; export function BoundaryConditionsPanel() { @@ -33,6 +34,7 @@ export function BoundaryConditionsPanel() { + )} diff --git a/web/src/components/panel/CouplingSection.tsx b/web/src/components/panel/CouplingSection.tsx new file mode 100644 index 0000000..4e46cb8 --- /dev/null +++ b/web/src/components/panel/CouplingSection.tsx @@ -0,0 +1,344 @@ +// SPDX-FileCopyrightText: 2026 Michael Kofler +// SPDX-License-Identifier: AGPL-3.0-or-later + +import { useEffect, useMemo, useState } from "react"; +import { useModelStore, referencePointOptions } from "../../store/modelStore"; +import type { + CouplingGroup, + CouplingKind, + ReferencePointOption, +} from "../../store/modelStore"; +import { usePickedFaces } from "../../hooks/usePickedFaces"; +import { DOF_LABELS, toFaceEntries } from "./bcFormUtils"; +import { + CouplingKindSelect, + DofCheckboxes, + PickedFaceList, + PickGeometryToggle, + ReferencePointInputs, +} from "./BcLoadFormControls"; +import { CouplingValueForm } from "./GroupValueForms"; +import { GroupCard } from "./GroupCard"; +import { fmt } from "../../lib/modelDisplay"; +import styles from "./LeftPanel.module.css"; + +// One-line summary of a coupling: how it ties, and — for a kinematic one — the +// DOFs it ties. A distributing coupling always ties all six, so listing them +// would say nothing. +export function couplingGroupMeta(group: CouplingGroup): string { + const nNodes = new Set(group.faces.flatMap((face) => face.nodeIds)).size; + const kind = group.kind === "kinematic" ? "kinematic" : "distributing"; + const dofs = + group.kind === "kinematic" + ? ` · ${group.dofs.map((dof) => DOF_LABELS[dof]).join(", ")}` + : ""; + return `${kind}${dofs} · ${nNodes} node${nNodes === 1 ? "" : "s"} → (${group.point + .map((c) => fmt(c)) + .join(", ")})`; +} + +// Parse the reference point's coordinate boxes. A point at the origin is +// perfectly valid, so only non-finite input is rejected — never coerced to 0, +// which would silently move the coupling somewhere the user did not ask for. +function parsePoint( + coords: [string, string, string], + onError: (msg: string) => void, +): [number, number, number] | null { + const parsed = coords.map((c) => parseFloat(c)); + if (parsed.some((c) => !isFinite(c))) { + onError("Each reference point coordinate must be a finite number"); + return null; + } + return [parsed[0], parsed[1], parsed[2]]; +} + +// Surface-to-point coupling section: pick the surface, choose how it ties to its +// reference point and where that point sits, then apply. The same shape as a BC +// or a load, because a coupling is model data just like they are. +export function CouplingSection({ + onError, +}: { + onError(msg: string | null): void; +}) { + const nodes = useModelStore((s) => s.nodes); + const elements = useModelStore((s) => s.elements); + const couplingGroups = useModelStore((s) => s.couplingGroups); + const pickMode = useModelStore((s) => s.pickMode); + const pickGeometry = useModelStore((s) => s.pickGeometry); + const setPickGeometry = useModelStore((s) => s.setPickGeometry); + // Edge picking only means something on a shell: a closed solid boundary has + // no boundary polyline to walk (extractBoundaryEdges finds none). + const hasShells = useModelStore((s) => + s.elements.some((el) => el.type === "CTRIA3"), + ); + const createCouplingGroup = useModelStore((s) => s.createCouplingGroup); + const addFaceToCouplingGroup = useModelStore((s) => s.addFaceToCouplingGroup); + const updateCouplingGroup = useModelStore((s) => s.updateCouplingGroup); + const removeFaceFromCouplingGroup = useModelStore( + (s) => s.removeFaceFromCouplingGroup, + ); + const deleteCouplingGroup = useModelStore((s) => s.deleteCouplingGroup); + const setCouplingDraft = useModelStore((s) => s.setCouplingDraft); + const { + pickTargetGroupId, + setPickMode, + allPickedFaces, + removePickedFace, + endPick, + startPickForGroup, + } = usePickedFaces(onError); + + const [editingId, setEditingId] = useState(null); + const [kind, setKind] = useState("kinematic"); + const [checkedDofs, setCheckedDofs] = useState([ + true, + true, + true, + true, + true, + true, + ]); + const [coords, setCoords] = useState<[string, string, string]>([ + "0", + "0", + "0", + ]); + // Which derived position the coordinates were last filled from, so the boxes + // follow the pick until the user types their own numbers. + const [placedLabel, setPlacedLabel] = useState(null); + + const targetGroup = + pickTargetGroupId !== null + ? (couplingGroups.find((group) => group.id === pickTargetGroupId) ?? null) + : null; + + // Positions derived from the surface picked so far — its centre, plus the two + // ends of its axis when it is a cylinder (KOF-208). + const options = useMemo( + () => referencePointOptions(allPickedFaces, nodes, elements), + [allPickedFaces, nodes, elements], + ); + // Default the point to the surface centre as soon as one is picked, and keep + // following it while more faces are added — until the user picks another + // position or types coordinates, at which point the boxes are theirs. + const centre = options.length > 0 ? options[0] : null; + const centreKey = centre ? centre.point.join(",") : ""; + const [followedCentre, setFollowedCentre] = useState(""); + if (centre && placedLabel === null && centreKey !== followedCentre) { + setFollowedCentre(centreKey); + setCoords([ + String(centre.point[0]), + String(centre.point[1]), + String(centre.point[2]), + ]); + } + + // Mirror the coupling being built into the store so the viewport can draw it: + // the point where it currently sits, and the nodes it would grip. A reference + // point is a position in space with nothing in the mesh to anchor it, so + // without this the coordinate boxes are the only feedback there is until the + // coupling is already applied. + const draftKey = `${pickMode}|${targetGroup?.id ?? ""}|${coords.join(",")}|${allPickedFaces + .map((face) => face.nodeIds.length) + .join(",")}`; + useEffect(() => { + // Before a surface is picked there is no coupling yet — the coordinate + // boxes still hold their initial zeros, and previewing those would put a + // marker at the origin that stands for nothing. (Adding faces to an existing + // coupling is not placing a point either; that group's own spider already + // shows where its point is.) + if (pickMode !== "coupling" || targetGroup || allPickedFaces.length === 0) { + setCouplingDraft(null); + return; + } + const parsed = coords.map((coord) => parseFloat(coord)); + if (parsed.some((coord) => !isFinite(coord))) { + // Half-typed coordinates are not a position — drop the preview rather than + // park the marker at a number the user did not mean. + setCouplingDraft(null); + return; + } + setCouplingDraft({ + point: [parsed[0], parsed[1], parsed[2]], + nodeIds: allPickedFaces.flatMap((face) => face.nodeIds), + }); + // Keyed on draftKey alone, deliberately: `coords` and `allPickedFaces` are + // rebuilt on every render, and re-running on those would set a fresh draft + // object each time, which re-renders, which re-runs — a loop. draftKey + // changes exactly when the preview should. + }, [draftKey, setCouplingDraft]); + + // The preview belongs to this form; it must not survive it. + useEffect(() => () => setCouplingDraft(null), [setCouplingDraft]); + + function placeAt(option: ReferencePointOption) { + setPlacedLabel(option.label); + setCoords([ + String(option.point[0]), + String(option.point[1]), + String(option.point[2]), + ]); + } + + function startCouplingPick(groupId: number | null) { + setPlacedLabel(null); + setFollowedCentre(""); + if (groupId === null) setPickMode("coupling", null); + else startPickForGroup("coupling", groupId); + } + + function applyCoupling() { + if (allPickedFaces.length === 0) { + onError("A coupling needs at least one picked face"); + return; + } + if (targetGroup) { + for (const faceEntry of toFaceEntries( + allPickedFaces, + targetGroup.faces.length, + pickGeometry, + )) + addFaceToCouplingGroup(targetGroup.id, faceEntry); + endPick(); + return; + } + const point = parsePoint(coords, onError); + if (point === null) return; + const dofs = checkedDofs + .map((checked, i) => (checked ? i : -1)) + .filter((i) => i >= 0); + if (kind === "kinematic" && dofs.length === 0) { + onError("A kinematic coupling must tie at least one DOF"); + return; + } + createCouplingGroup( + toFaceEntries(allPickedFaces, 0, pickGeometry), + point, + kind, + dofs, + ); + endPick(); + } + + return ( + <> +
+ Couplings +
+ + {pickMode !== "coupling" && ( + + )} + + {pickMode === "coupling" && ( +
+
+ + {targetGroup ? `Add face to ${targetGroup.name}` : "New Coupling"} + + +
+ + {/* A coupling can grip a LINE as well as a surface — the rim of a + shell, a stiffener edge — which on a flat sheet is the only way to + select it at all, since the whole sheet is one face-pick region. */} + + + + + {!targetGroup && allPickedFaces.length > 0 && ( + <> + + {/* The DOF mask only exists for a kinematic coupling — a + distributing one ties all six of its reference point's DOFs by + construction, so offering checkboxes would promise a control + the solver does not have. */} + {kind === "kinematic" && ( + + setCheckedDofs((prev) => + prev.map((checked, i) => + i === index ? !checked : checked, + ), + ) + } + /> + )} + { + setPlacedLabel("custom"); + setCoords((prev) => { + const next = [...prev] as [string, string, string]; + next[index] = value; + return next; + }); + }} + onPickOption={placeAt} + /> + + )} + + +
+ )} + + {couplingGroups.map((group) => ( + { + updateCouplingGroup(group.id, nextKind, nextDofs, nextPoint); + setEditingId(null); + }} + onCancel={() => setEditingId(null)} + /> + ) + } + onStartPick={() => startCouplingPick(group.id)} + onToggleEdit={() => + setEditingId(editingId === group.id ? null : group.id) + } + onDelete={() => deleteCouplingGroup(group.id)} + onRemoveFace={(faceId) => + removeFaceFromCouplingGroup(group.id, faceId) + } + /> + ))} + + ); +} diff --git a/web/src/components/panel/GroupValueForms.tsx b/web/src/components/panel/GroupValueForms.tsx index b5d2070..9f3cfa6 100644 --- a/web/src/components/panel/GroupValueForms.tsx +++ b/web/src/components/panel/GroupValueForms.tsx @@ -1,13 +1,16 @@ // SPDX-FileCopyrightText: 2026 Michael Kofler // SPDX-License-Identifier: AGPL-3.0-or-later -import { useState } from "react"; +import { useEffect, useMemo, useState } from "react"; import { loadKind, loadComponents, + referencePointOptions, useModelStore, } from "../../store/modelStore"; import type { + CouplingGroup, + CouplingKind, LoadKind, NamedBcGroup, NamedLoadGroup, @@ -21,10 +24,12 @@ import { parseTieDistance, } from "./bcFormUtils"; import { + CouplingKindSelect, DofCheckboxes, LoadKindSelect, LoadVectorInputs, PressureInput, + ReferencePointInputs, TieExtentInputs, } from "./BcLoadFormControls"; import styles from "./LeftPanel.module.css"; @@ -236,3 +241,128 @@ export function LoadValueForm({ ); } + +// Inline editor for a coupling's kind, tied DOFs and reference point position — +// opened by the coupling's ✎ button. The gripped surface is edited through the +// face rows, exactly as for a BC or a load. +export function CouplingValueForm({ + group, + onSave, + onCancel, +}: { + group: CouplingGroup; + onSave( + kind: CouplingKind, + dofs: number[], + point: [number, number, number], + ): void; + onCancel(): void; +}) { + const nodes = useModelStore((s) => s.nodes); + const elements = useModelStore((s) => s.elements); + const setCouplingDraft = useModelStore((s) => s.setCouplingDraft); + const [kind, setKind] = useState(group.kind); + const [checkedDofs, setCheckedDofs] = useState( + DOF_LABELS.map((_, i) => group.dofs.includes(i)), + ); + const [coords, setCoords] = useState<[string, string, string]>([ + String(group.point[0]), + String(group.point[1]), + String(group.point[2]), + ]); + const [error, setError] = useState(null); + + // Re-derived from the coupling's own surface, so a point can be re-centred + // after the surface was extended or trimmed. + const options = useMemo( + () => referencePointOptions(group.faces, nodes, elements), + [group.faces, nodes, elements], + ); + + // Preview the edited position in the viewport, alongside the coupling's + // committed spider — so the point being moved is visible next to where it + // currently is, which is the comparison the edit is about. + const draftKey = coords.join(","); + useEffect(() => { + const parsed = coords.map((coord) => parseFloat(coord)); + if (parsed.some((coord) => !isFinite(coord))) { + setCouplingDraft(null); + return; + } + setCouplingDraft({ + point: [parsed[0], parsed[1], parsed[2]], + nodeIds: group.faces.flatMap((face) => face.nodeIds), + }); + // Keyed on draftKey rather than `coords`, which is a fresh array on every + // render: re-running on it would set a new draft object each time, which + // re-renders, which re-runs. + }, [draftKey, group.faces, setCouplingDraft]); + + // The preview belongs to this form; closing it must take the preview too. + useEffect(() => () => setCouplingDraft(null), [setCouplingDraft]); + + function handleSave() { + const parsed = coords.map((c) => parseFloat(c)); + if (parsed.some((c) => !isFinite(c))) { + setError("Each reference point coordinate must be a finite number"); + return; + } + const dofs = checkedDofs + .map((checked, i) => (checked ? i : -1)) + .filter((i) => i >= 0); + if (kind === "kinematic" && dofs.length === 0) { + setError("A kinematic coupling must tie at least one DOF"); + return; + } + onSave(kind, dofs, [parsed[0], parsed[1], parsed[2]]); + } + + return ( +
+ {error && ( +
+ {error} + +
+ )} + + {kind === "kinematic" && ( + + setCheckedDofs((prev) => + prev.map((checked, i) => (i === index ? !checked : checked)), + ) + } + /> + )} + + setCoords((prev) => { + const next = [...prev] as [string, string, string]; + next[index] = value; + return next; + }) + } + onPickOption={(option) => + setCoords([ + String(option.point[0]), + String(option.point[1]), + String(option.point[2]), + ]) + } + /> +
+ + +
+
+ ); +} diff --git a/web/src/components/panel/LeftPanel.module.css b/web/src/components/panel/LeftPanel.module.css index de08713..dfa5d14 100644 --- a/web/src/components/panel/LeftPanel.module.css +++ b/web/src/components/panel/LeftPanel.module.css @@ -841,6 +841,16 @@ flex-shrink: 0; } +/* Surface-to-point couplings — the emerald of the coupling spider in the + viewport, so a coupling reads the same in the panel and in 3D. */ +.couplingDot { + width: 8px; + height: 8px; + border-radius: 50%; + background: #059669; + flex-shrink: 0; +} + /* Header of one side (A / B) of a connection's face list. */ .tieSideRow { display: flex; diff --git a/web/src/components/panel/LoadSection.tsx b/web/src/components/panel/LoadSection.tsx index e022c9f..8f946c2 100644 --- a/web/src/components/panel/LoadSection.tsx +++ b/web/src/components/panel/LoadSection.tsx @@ -39,6 +39,15 @@ export function LoadSection({ const hasShells = useModelStore((s) => s.elements.some((el) => el.type === "CTRIA3"), ); + // A nodal MOMENT only reaches the solver on the pure-shell path + // (shellPointLoads carries DOFs 3..5). The solid and coupled assemblers give + // an ordinary mesh node three translations only and drop a rotational load, + // so on those a couple has to go through a coupling's reference point. + const isPureShell = useModelStore( + (s) => + s.elements.length > 0 && s.elements.every((el) => el.type === "CTRIA3"), + ); + const couplingGroups = useModelStore((s) => s.couplingGroups); const createLoadGroup = useModelStore((s) => s.createLoadGroup); const addFaceToLoadGroup = useModelStore((s) => s.addFaceToLoadGroup); const removeFaceFromLoadGroup = useModelStore( @@ -72,6 +81,17 @@ export function LoadSection({ "1000", ]); const [pressureVal, setPressureVal] = useState("10"); + // Whether every picked selection is a coupling's REFERENCE POINT — clicked in + // the viewport like anything else. Such a point carries rotational DOFs, so + // it is the one place a MOMENT can act on a solid model: the couple is + // applied to the point and the coupling spreads it over the gripped surface + // (rebuildLoads / coupledLoads). + const refPoints = new Set(couplingGroups.map((c) => c.refNodeId)); + const onlyReferencePoints = + allPickedFaces.length > 0 && + allPickedFaces.every((face) => + face.nodeIds.every((id) => refPoints.has(id)), + ); const targetLoadGroup = pickTargetGroupId !== null @@ -84,17 +104,43 @@ export function LoadSection({ allPickedFaces, // eslint-disable-next-line kofem/no-silent-fallback -- numbering offset for the new entries; a pick with no target group starts a fresh group, which has 0 faces targetLoadGroup?.faces.length ?? 0, - pickGeometry === "edge" ? "Edge" : "Face", + pickGeometry, ); if (targetLoadGroup) { for (const faceEntry of faceEntries) { addFaceToLoadGroup(targetLoadGroup.id, faceEntry); } } else if (loadKindSel === "pressure") { + // A pressure is force per unit area, and a single node — a mesh node or a + // reference point alike — has none. Refuse it here rather than create a + // group whose selection contributes nothing to the solve. + if (pickGeometry === "point") { + onError( + "A pressure cannot act on a single node — it has no area. Apply a force there, or pick a face.", + ); + return; + } const pressure = parsePressure(pressureVal, onError); if (pressure === null) return; createLoadGroup(faceEntries, 0, pressure, "pressure"); } else { + // A moment needs a rotational DOF to act on. A coupling's reference point + // has one on every solve path; an ordinary mesh node has one only on the + // pure-shell path (shellPointLoads carries DOFs 3..5), while the solid and + // coupled assemblers would drop it. + if ( + loadKindSel === "moment" && + pickGeometry === "point" && + !onlyReferencePoints && + !isPureShell + ) { + onError( + "A moment cannot act on a single node of a solid mesh — its nodes carry no " + + "rotational DOF, so the couple would be dropped. Couple the surface to a " + + "reference point and apply the moment there.", + ); + return; + } const components = parseLoadVector( loadKindSel === "moment" ? momentVec : forceVec, loadKindSel, @@ -106,6 +152,31 @@ export function LoadSection({ endPick(); } + // How the load about to be created will actually reach the solver. The three + // routes are genuinely different physics — an integrated traction, a + // concentrated nodal load, and a couple carried by a coupling — so which one + // a selection lands on is worth saying before it is applied. + function applicationNote(): string { + if (loadKindSel === "pressure") + return "applied as p·n̂ over each face (work-equivalent)"; + if (loadKindSel === "force") { + if (onlyReferencePoints) + return "applied at the reference point, and spread over its coupled surface"; + if (pickGeometry === "point") + return "applied at the picked node as a concentrated force"; + if (pickGeometry === "edge") + return "applied as a work-equivalent line load along the edge"; + return "applied as a work-equivalent surface traction"; + } + if (onlyReferencePoints) + return "applied to the reference point as a couple, and spread over its coupled surface"; + if (pickGeometry === "point") + return isPureShell + ? "applied at the picked node as a couple" + : "a solid mesh node carries no rotation — use a coupling's reference point"; + return "distributed as equivalent nodal forces"; + } + function updateVecComponent(index: number, value: string) { const setVec = loadKindSel === "moment" ? setMomentVec : setForceVec; setVec((prev) => { @@ -143,12 +214,15 @@ export function LoadSection({ - {hasShells && ( - - )} + {/* A point load is applied at the node itself rather than integrated + over a surface, which is the only way to state a concentrated + force — a lug pin, a bolt reaction — on a mesh that has no element + face to spread it over. */} + )} -
- {loadKindSel === "pressure" - ? "applied as p·n̂ over each face (work-equivalent)" - : loadKindSel === "force" - ? pickGeometry === "edge" - ? "applied as a work-equivalent line load along the edge" - : "applied as a work-equivalent surface traction" - : "distributed as equivalent nodal forces"} -
+
{applicationNote()}
diff --git a/web/src/components/panel/bcFormUtils.ts b/web/src/components/panel/bcFormUtils.ts index 6b5e813..7f27026 100644 --- a/web/src/components/panel/bcFormUtils.ts +++ b/web/src/components/panel/bcFormUtils.ts @@ -8,6 +8,8 @@ import type { NamedLoadGroup, } from "../../store/modelStore"; import { fmt } from "../../lib/modelDisplay"; +import { SELECTION_NOUN } from "../../lib/facePick"; +import type { PickGeometry } from "../../lib/facePick"; export const DOF_LABELS = ["Ux", "Uy", "Uz", "Rx", "Ry", "Rz"]; export const FORCE_LABELS = ["Fx", "Fy", "Fz"]; @@ -32,14 +34,21 @@ export function faceKey(face: FaceSelection): string { return `${face.nodeIds.length}-${face.nodeIds[0]}-${face.nodeIds[face.nodeIds.length - 1]}`; } +// Committed entries for a set of picked selections. The picked GEOMETRY rides +// along on each entry: only a "point" changes how the solver applies the group +// (a single node spans no element face, so a load on it is applied at the node +// rather than integrated as a traction — see rebuildLoads), and the panel must +// state which it was rather than let a downstream builder guess from the node +// count. export function toFaceEntries( faces: FaceSelection[], existingCount: number, - noun: "Face" | "Edge" = "Face", + geometry: PickGeometry = "face", ) { return faces.map((face, i) => ({ - label: `${noun} ${existingCount + i + 1}`, + label: `${SELECTION_NOUN[geometry]} ${existingCount + i + 1}`, nodeIds: face.nodeIds, + geometry, })); } diff --git a/web/src/components/statusbar/StatusBar.tsx b/web/src/components/statusbar/StatusBar.tsx index 426a79f..6ff1356 100644 --- a/web/src/components/statusbar/StatusBar.tsx +++ b/web/src/components/statusbar/StatusBar.tsx @@ -2,16 +2,18 @@ // SPDX-License-Identifier: AGPL-3.0-or-later import { useModelStore, loadKind } from "../../store/modelStore"; +import type { PickMode } from "../../store/modelStore"; import { APP_VERSION } from "../../lib/version"; import styles from "./StatusBar.module.css"; type ViewRepr = "geometry" | "surface" | "volume" | "wireframe"; // What the face being picked is for, shown next to the pick indicator. -const PICK_PURPOSE: Record<"bc" | "load" | "tie", string> = { +const PICK_PURPOSE: Record = { bc: "fixed displacement", load: "apply load", tie: "tie connection", + coupling: "surface-to-point coupling", }; const REPR_BUTTONS: { diff --git a/web/src/components/viewport/BoundaryConditionLayer.tsx b/web/src/components/viewport/BoundaryConditionLayer.tsx index 7d48269..1deac0a 100644 --- a/web/src/components/viewport/BoundaryConditionLayer.tsx +++ b/web/src/components/viewport/BoundaryConditionLayer.tsx @@ -1,10 +1,11 @@ // SPDX-FileCopyrightText: 2026 Michael Kofler // SPDX-License-Identifier: AGPL-3.0-or-later -// Boundary condition, load and tie visualisation: committed BC/load/tie face -// highlights, the in-progress pick-session highlights (pending + selected -// faces), the fixed-support marker, resultant / per-node load arrows, and the -// links between the node pairs a tie connection welds. +// Boundary condition, load, tie and coupling visualisation: committed +// BC/load/tie/coupling face highlights, the in-progress pick-session highlights +// (pending + selected faces), the fixed-support marker, resultant / per-node +// load arrows, the links between the node pairs a tie connection welds, and the +// spider of a surface-to-point coupling. import { useMemo } from "react"; import * as THREE from "three"; @@ -14,6 +15,7 @@ import { loadComponents, } from "../../store/modelStore"; import { useTiePairs } from "../../hooks/useTiePairs"; +import { useReferencePointPick } from "./useFacePick"; import { buildFacePositions } from "./useMeshTopology"; import type { MeshTopology } from "./useMeshTopology"; @@ -22,6 +24,17 @@ import type { MeshTopology } from "./useMeshTopology"; const TIE_COLOR_A = "#7c3aed"; const TIE_COLOR_B = "#0891b2"; +// Couplings: emerald, matching the panel's coupling dot. A coupling still being +// placed gets its own blue rather than the pick session's orange or the +// committed emerald: it sits ON the orange picked surface, and when a point is +// being MOVED it has to be told apart from the committed spider still drawn at +// the old position — which is the whole comparison the edit is about. +const COUPLING_COLOR = "#059669"; +const COUPLING_DRAFT_COLOR = "#2563eb"; + +// The orange every live pick session uses for what is currently selected. +const SELECTION_COLOR = "#e05533"; + interface BoundaryConditionLayerProps { topology: MeshTopology; showResult: boolean; @@ -39,10 +52,24 @@ export function BoundaryConditionLayer({ const bcGroups = useModelStore((s) => s.bcGroups); const loadGroups = useModelStore((s) => s.loadGroups); const tieGroups = useModelStore((s) => s.tieGroups); + const couplingGroups = useModelStore((s) => s.couplingGroups); + const couplingDraft = useModelStore((s) => s.couplingDraft); const pickMode = useModelStore((s) => s.pickMode); const tieDraft = useModelStore((s) => s.tieDraft); const loadDisplay = useModelStore((s) => s.loadDisplay); const { pairs: tiePairs } = useTiePairs(); + // Set only while a BC or load is being picked in point mode — see + // useReferencePointPick. When it is, the markers are click targets and are + // drawn larger to say so. + const pickReferencePoint = useReferencePointPick(); + // Nodes in the live pick session, so a selected reference point is shown + // selected rather than the user having to read it off the panel. + const pickedNodeIds = useMemo(() => { + const ids = new Set(); + for (const face of pendingFaces) for (const id of face.nodeIds) ids.add(id); + if (selectedFace) for (const id of selectedFace.nodeIds) ids.add(id); + return ids; + }, [pendingFaces, selectedFace]); const { nodeMap, @@ -176,6 +203,80 @@ export function BoundaryConditionLayer({ return { positions: positions.subarray(0, 6 * written), markers }; }, [tiePairs, nodeMap]); + // Coupling surface highlights — the surface each coupling grips. + const couplingFaceHighlights = useMemo(() => { + if (!boundaryMeshTopo) return null; + const highlights: { key: string; positions: Float32Array }[] = []; + for (const group of couplingGroups) + group.faces.forEach((face, i) => { + const positions = buildFacePositions( + face.nodeIds, + boundaryMeshTopo.triangles, + nodeMap, + ); + if (positions) + highlights.push({ key: `coupling-${group.id}-${i}`, positions }); + }); + return highlights; + }, [couplingGroups, boundaryMeshTopo, nodeMap]); + + // Coupling spiders — a line from every gripped node to its reference point, + // plus a marker at the point. This is the constraint's actual reach: which + // nodes the point governs (kinematic) or averages (distributing), rather than + // just the surface it was declared on. The reference point is a node of the + // model, so its position comes from the same nodeMap as everything else. + const couplingSpiders = useMemo(() => { + const spiders: { + id: number; + refNodeId: number; + positions: Float32Array; + point: [number, number, number]; + }[] = []; + for (const group of couplingGroups) { + const ref = nodeMap.get(group.refNodeId); + if (!ref) continue; + const gripped = new Set(group.faces.flatMap((face) => face.nodeIds)); + gripped.delete(group.refNodeId); + const positions = new Float32Array(6 * gripped.size); + let written = 0; + for (const nodeId of gripped) { + const node = nodeMap.get(nodeId); + if (!node) continue; + positions.set( + [ref.n.x, ref.n.y, ref.n.z, node.n.x, node.n.y, node.n.z], + 6 * written, + ); + written++; + } + spiders.push({ + id: group.id, + refNodeId: group.refNodeId, + positions: positions.subarray(0, 6 * written), + point: [ref.n.x, ref.n.y, ref.n.z], + }); + } + return spiders; + }, [couplingGroups, nodeMap]); + + // The coupling being placed: the same spider, at the position the form + // currently holds. Drawn from `couplingDraft` rather than re-derived here, so + // the marker is exactly where the coordinates say and cannot drift from them. + const draftSpider = useMemo(() => { + if (!couplingDraft) return null; + const positions = new Float32Array(6 * couplingDraft.nodeIds.length); + let written = 0; + for (const nodeId of new Set(couplingDraft.nodeIds)) { + const node = nodeMap.get(nodeId); + if (!node) continue; + positions.set( + [...couplingDraft.point, node.n.x, node.n.y, node.n.z], + 6 * written, + ); + written++; + } + return { positions: positions.subarray(0, 6 * written) }; + }, [couplingDraft, nodeMap]); + // BC markers — one small triangular cone per constrained node (apex at the // node, base outward), replacing the former single centroid marker. Each // marker is oriented along the outward normal of the constrained surface, @@ -583,6 +684,107 @@ export function BoundaryConditionLayer({ )} + {/* Coupling surface highlights — the surface each coupling grips */} + {!showResult && + couplingFaceHighlights?.map((highlight) => ( + + + + + + + ))} + + {/* Coupling spiders — a line from the reference point to every node it + couples, and a sphere at the point itself */} + {!showResult && + couplingSpiders.map((spider) => ( + + {spider.positions.length > 0 && ( + + + + + + + )} + {/* `transparent` with full opacity, not a translucency: three.js + draws the opaque pass before the transparent one, so an opaque + marker is painted UNDER the surface highlights however high its + renderOrder. Joining the transparent pass is what lets the order + put the point on top of the surface it belongs to. */} + { + event.stopPropagation(); + pickReferencePoint(spider.refNodeId); + }) + } + > + + + + + ))} + + {/* The coupling being placed — its point and the nodes it would grip, + so the position can be judged against the model before it is applied */} + {!showResult && couplingDraft && draftSpider && ( + + {draftSpider.positions.length > 0 && ( + + + + + + + )} + + + + + + )} + {/* Pending faces — accumulated via shift-click, same colour as selection but slightly dimmer */} {pendingFacePositions && ( diff --git a/web/src/components/viewport/useFacePick.ts b/web/src/components/viewport/useFacePick.ts index 1c387a3..5741c56 100644 --- a/web/src/components/viewport/useFacePick.ts +++ b/web/src/components/viewport/useFacePick.ts @@ -6,7 +6,10 @@ import { useModelStore } from "../../store/modelStore"; import { pickFaceNodeIds, pickEdgeNodeIds, + pickPointNodeId, + nearestReferencePoint, toggleFaceSelection, + SELECTION_NOUN, } from "../../lib/facePick"; import type { BoundaryMeshTopo, Vec3 } from "../../lib/facePick"; @@ -19,6 +22,7 @@ import type { BoundaryMeshTopo, Vec3 } from "../../lib/facePick"; // otherwise a BFS flood-fill with normal-angle thresholds (parametric box mesh // or .inp import). Edge mode selects the boundary polyline near the click — the // only way to grab the rim of a flat shell, whose whole sheet is one region. +// Point mode selects the single nearest node of the clicked facet. export function useFacePick( boundaryMeshTopo: BoundaryMeshTopo | null, getPos: (id: number) => Vec3, @@ -29,6 +33,7 @@ export function useFacePick( const pendingFaces = useModelStore((s) => s.pendingFaces); const setSelectedFace = useModelStore((s) => s.setSelectedFace); const setPendingFaces = useModelStore((s) => s.setPendingFaces); + const couplingGroups = useModelStore((s) => s.couplingGroups); function handleFacePick(e: ThreeEvent) { if (!pickMode || e.faceIndex == null || !boundaryMeshTopo) return; @@ -37,18 +42,38 @@ export function useFacePick( const startIdx = e.faceIndex; if (startIdx >= boundaryMeshTopo.triangles.length) return; - const faceNodeIds = + const clickPoint: Vec3 = [e.point.x, e.point.y, e.point.z]; + const picked = pickGeometry === "edge" - ? [ - ...pickEdgeNodeIds( - [e.point.x, e.point.y, e.point.z], - startIdx, - boundaryMeshTopo, - getPos, - ), - ] - : [...pickFaceNodeIds(startIdx, boundaryMeshTopo)]; + ? pickEdgeNodeIds(clickPoint, startIdx, boundaryMeshTopo, getPos) + : pickGeometry === "point" + ? pickPointNodeId(clickPoint, startIdx, boundaryMeshTopo, getPos) + : pickFaceNodeIds(startIdx, boundaryMeshTopo); + let faceNodeIds = [...picked]; if (faceNodeIds.length === 0) return; + + // In point mode a coupling's REFERENCE POINT competes with the mesh node. + // It belongs to no surface, so the ray can only ever report the mesh behind + // it — which on a marker drawn over the face it couples is always the mesh. + // Whether the point or the node was meant is decided by which is nearer to + // where the click landed, the same rule pickPointNodeId uses to choose + // between the corners of a facet. + if (pickGeometry === "point" && couplingGroups.length > 0) { + const nearestNode = getPos(faceNodeIds[0]); + const nodeDistSq = + (clickPoint[0] - nearestNode[0]) ** 2 + + (clickPoint[1] - nearestNode[1]) ** 2 + + (clickPoint[2] - nearestNode[2]) ** 2; + const nearestRef = nearestReferencePoint( + clickPoint, + couplingGroups.map((group) => ({ + nodeId: group.refNodeId, + point: group.point, + })), + ); + if (nearestRef && nearestRef.distSq < nodeDistSq) + faceNodeIds = [nearestRef.nodeId]; + } // The normal decides which axis/extreme the picked face is snapped to, and // that selection becomes solver input. Substituting +Y for a missing normal // would silently snap the pick to the wrong face — drop the pick instead. @@ -84,7 +109,7 @@ export function useFacePick( isMax, label: "", }, - pickGeometry === "edge" ? "Edge" : "Face", + SELECTION_NOUN[pickGeometry], ); // Re-split into pending faces + the active (last) selection. @@ -94,3 +119,44 @@ export function useFacePick( return pickMode ? handleFacePick : undefined; } + +// Picking handler for a coupling's REFERENCE POINT marker — the second of the +// two routes by which a reference point is selected. +// +// The first is the distance test inside handleFacePick above, which covers the +// usual case: the marker is drawn on or near the surface it couples, so the ray +// reports that surface and the point is chosen by being nearer to the click. +// This one covers the case that test cannot reach — a point placed clear of the +// model, where the ray hits nothing at all and no surface event is raised. +// +// Only in POINT mode: in face or edge mode the marker would sit in front of the +// surface being picked and swallow clicks meant for it. Returns undefined +// otherwise, so the marker carries no onClick at all and stays inert. +export function useReferencePointPick(): + ((nodeId: number) => void) | undefined { + const pickMode = useModelStore((s) => s.pickMode); + const pickGeometry = useModelStore((s) => s.pickGeometry); + const selectedFace = useModelStore((s) => s.selectedFace); + const pendingFaces = useModelStore((s) => s.pendingFaces); + const setSelectedFace = useModelStore((s) => s.setSelectedFace); + const setPendingFaces = useModelStore((s) => s.setPendingFaces); + + function pickReferencePoint(nodeId: number) { + const current = selectedFace + ? [...pendingFaces, selectedFace] + : pendingFaces; + // A reference point has no surface and so no normal; the axis/isMax fields + // exist to snap a FACE pick to an extreme and mean nothing here. + const next = toggleFaceSelection( + current, + { nodeIds: [nodeId], axis: "X", isMax: true, label: "" }, + SELECTION_NOUN.point, + ); + setPendingFaces(next.slice(0, -1)); + setSelectedFace(next.length > 0 ? next[next.length - 1] : null); + } + + return (pickMode === "bc" || pickMode === "load") && pickGeometry === "point" + ? pickReferencePoint + : undefined; +} diff --git a/web/src/hooks/usePickedFaces.ts b/web/src/hooks/usePickedFaces.ts index a5f13b4..94ab24c 100644 --- a/web/src/hooks/usePickedFaces.ts +++ b/web/src/hooks/usePickedFaces.ts @@ -2,6 +2,7 @@ // SPDX-License-Identifier: AGPL-3.0-or-later import { useModelStore } from "../store/modelStore"; +import type { PickMode } from "../store/modelStore"; // Shared face-picking session state: the viewport writes clicked faces into // the store (selectedFace + shift-click pendingFaces); both pick panels read @@ -34,7 +35,7 @@ export function usePickedFaces(onError: (msg: string | null) => void) { setPendingFaces([]); } - function startPickForGroup(kind: "bc" | "load" | "tie", groupId: number) { + function startPickForGroup(kind: PickMode, groupId: number) { setPickMode(kind, groupId); setSelectedFace(null); setPendingFaces([]); diff --git a/web/src/hooks/useSolver.ts b/web/src/hooks/useSolver.ts index 173845b..bfacd9a 100644 --- a/web/src/hooks/useSolver.ts +++ b/web/src/hooks/useSolver.ts @@ -25,6 +25,7 @@ export function useSolver() { const setMode = useModelStore((s) => s.setMode); const elementOrder = useModelStore((s) => s.elementOrder); const tieGroups = useModelStore((s) => s.tieGroups); + const couplingGroups = useModelStore((s) => s.couplingGroups); // Surface mesh + CAD face ids drive auto-shell idealisation of thin bodies // in the worker's coupled solve. const surfaceTriangles = useModelStore((s) => s.surfaceTriangles); @@ -67,6 +68,7 @@ export function useSolver() { surfaceLoads, elementOrder, tieGroups, + couplings: couplingGroups, surfaceTriangles, surfaceFaceIds, }, diff --git a/web/src/lib/analysisFile.ts b/web/src/lib/analysisFile.ts index 629f00d..69ce781 100644 --- a/web/src/lib/analysisFile.ts +++ b/web/src/lib/analysisFile.ts @@ -3,6 +3,7 @@ import type { AppMode, + CouplingGroup, Element, ElementType, Material, @@ -72,9 +73,11 @@ export interface AnalysisState { bcGroups: NamedBcGroup[]; loadGroups: NamedLoadGroup[]; tieGroups: TieGroup[]; + couplingGroups: CouplingGroup[]; nextBcGroupId: number; nextLoadGroupId: number; nextTieGroupId: number; + nextCouplingGroupId: number; nextFaceEntryId: number; nextMatId: number; stepSurface: StepTessellation | null; @@ -104,9 +107,14 @@ interface KofemFieldDataV1 { // Tie (connector) conditions. Absent in files written before connections // became a named model object — such a file simply has no ties. tieGroups?: TieGroup[]; + // Surface-to-point couplings. Absent in files written before couplings + // existed — such a file simply has none. The reference points themselves need + // no separate entry: they are nodes, and travel in the VTU point list. + couplingGroups?: CouplingGroup[]; nextBcGroupId: number; nextLoadGroupId: number; nextTieGroupId?: number; + nextCouplingGroupId?: number; nextFaceEntryId: number; nextMatId: number; stepSurface: StepTessellation | null; @@ -242,9 +250,11 @@ export function serializeAnalysis(state: AnalysisState): string { bcGroups: state.bcGroups, loadGroups: state.loadGroups, tieGroups: state.tieGroups, + couplingGroups: state.couplingGroups, nextBcGroupId: state.nextBcGroupId, nextLoadGroupId: state.nextLoadGroupId, nextTieGroupId: state.nextTieGroupId, + nextCouplingGroupId: state.nextCouplingGroupId, nextFaceEntryId: state.nextFaceEntryId, nextMatId: state.nextMatId, stepSurface: state.stepSurface, @@ -403,6 +413,10 @@ function parseMetadata(xml: string): KofemFieldDataV1 { throw new Error( `Invalid analysis file: "tieGroups" must be an array, got ${typeof meta.tieGroups}`, ); + if (meta.couplingGroups !== undefined && !Array.isArray(meta.couplingGroups)) + throw new Error( + `Invalid analysis file: "couplingGroups" must be an array, got ${typeof meta.couplingGroups}`, + ); if (typeof meta.modelName !== "string") throw new Error('Invalid analysis file: "modelName" must be a string'); @@ -523,6 +537,10 @@ export function parseAnalysisFile(text: string): AnalysisState { const nextTieGroupId = meta.nextTieGroupId ?? tieGroups.reduce((highest, tie) => Math.max(highest, tie.id), 0) + 1; + const couplingGroups = meta.couplingGroups ?? []; + const nextCouplingGroupId = + meta.nextCouplingGroupId ?? + couplingGroups.reduce((highest, c) => Math.max(highest, c.id), 0) + 1; return { modelName: meta.modelName, @@ -535,9 +553,11 @@ export function parseAnalysisFile(text: string): AnalysisState { bcGroups: meta.bcGroups, loadGroups: meta.loadGroups, tieGroups, + couplingGroups, nextBcGroupId: meta.nextBcGroupId, nextLoadGroupId: meta.nextLoadGroupId, nextTieGroupId, + nextCouplingGroupId, nextFaceEntryId: meta.nextFaceEntryId, nextMatId: meta.nextMatId, stepSurface: meta.stepSurface, diff --git a/web/src/lib/coupling.ts b/web/src/lib/coupling.ts new file mode 100644 index 0000000..8d1955c --- /dev/null +++ b/web/src/lib/coupling.ts @@ -0,0 +1,464 @@ +// SPDX-FileCopyrightText: 2026 Michael Kofler +// SPDX-License-Identifier: AGPL-3.0-or-later + +// Surface-to-point coupling (KOF-208): a picked surface is idealised to a single +// REFERENCE POINT, so a bolt, a bearing or a screw connection is stated as the +// one point it acts through instead of being smeared over the nodes of a hole. +// +// The engine half landed first (engine/cpp/shell_core.cpp); this is the model +// side that feeds it. Two kinds, and the choice between them is a modelling +// decision, not a detail: +// +// distributing (RBE3) — the reference point is DEPENDENT: it follows the +// weighted average of the coupled nodes. It adds no +// stiffness, so the gripped surface stays flexible — +// but the point can only be LOADED, never fixed, +// driven, or coupled onward. +// kinematic (RBE2) — the coupled nodes are dependent and follow the point +// rigidly, u_i = u_R + θ_R × r_i. The point stays +// independent, so it can be fixed, loaded or tied to +// another point; the price is that the surface it grips +// becomes rigid. +// +// The DOF mask selects which of the reference point's six DOFs a KINEMATIC +// coupling ties. On a surface of solid nodes it does almost nothing (solid nodes +// have no rotational DOF to tie, and the point's rotation already reaches them +// through θ_R × r); it earns its keep on a point-to-point coupling, where +// x,y,z alone is a spherical joint and all six is a rigid link. + +// ── Model types ─────────────────────────────────────────────────────────────── + +export type CouplingKind = "distributing" | "kinematic"; + +export interface CouplingNode { + id: number; + x: number; + y: number; + z: number; +} + +export interface CouplingElement { + type: string; + nodeIds: number[]; +} + +// The structural minimum of a coupling the solver needs. The store's +// CouplingGroup (boundarySlice) satisfies it, and so does the solve payload. +export interface CouplingDefinition { + name: string; + kind: CouplingKind; + // Tied DOFs (0..2 translations, 3..5 rotations) of a kinematic coupling. + // Ignored by a distributing coupling, which always ties all six. + dofs: number[]; + // The reference point, as a node in the model's own node numbering. It is a + // real node (created with the coupling, removed with it), so a BC or a load + // reaches it through exactly the machinery every other node uses. + refNodeId: number; + // The coupled surface. + faces: { nodeIds: number[] }[]; +} + +// All six DOFs — what a coupling with no explicit mask ties. +export const ALL_DOFS = [0, 1, 2, 3, 4, 5]; + +// Engine `mpc` kind code (solve_coupled): 0 distributing RBE3, 2 kinematic RBE2. +// 1 is the shell-to-solid relaxed MPC, which is not a user-declared coupling. +export function couplingMpcCode(kind: CouplingKind): number { + return kind === "kinematic" ? 2 : 0; +} + +// Bit c of the engine's dof_mask selects DOF c of every coupled node. +export function couplingDofMask(kind: CouplingKind, dofs: number[]): number { + if (kind !== "kinematic") return 0x3f; // distributing ties all six + let mask = 0; + for (const dof of dofs) mask |= 1 << dof; + return mask; +} + +// The coupled nodes of one coupling, de-duplicated. A surface picked in two +// faces that overlap would otherwise offer the same node twice, and a kinematic +// coupling eliminating the same DOF twice is an error the engine refuses. +export function coupledNodeIds(coupling: CouplingDefinition): number[] { + const ids = new Set(); + for (const face of coupling.faces) + for (const nodeId of face.nodeIds) ids.add(nodeId); + ids.delete(coupling.refNodeId); + return [...ids]; +} + +// ── Reference point placement ───────────────────────────────────────────────── + +export interface ReferencePointOption { + // What the position means, for the panel's dropdown. + label: string; + point: [number, number, number]; +} + +type Vec3 = [number, number, number]; + +function sub(a: Vec3, b: Vec3): Vec3 { + return [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; +} + +function dot(a: Vec3, b: Vec3): number { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +function cross(a: Vec3, b: Vec3): Vec3 { + return [ + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ]; +} + +// Component-wise mean. Summed first and divided once: dividing each term by the +// count instead accumulates rounding, so the centre of a face whose coordinates +// average exactly comes back a few ulp off — visible as 1.9999999999999998 in a +// coordinate box the user is meant to read and re-type. +function centroidOf(points: Vec3[]): Vec3 { + const total: Vec3 = [0, 0, 0]; + for (const point of points) { + total[0] += point[0]; + total[1] += point[1]; + total[2] += point[2]; + } + return [ + total[0] / points.length, + total[1] / points.length, + total[2] / points.length, + ]; +} + +// Eigenvector of the smallest eigenvalue of a symmetric 3×3 matrix, by cyclic +// Jacobi rotations. The matrix is tiny and the sweep converges in a handful of +// passes, so this is exact enough to fit an axis and needs no linear-algebra +// dependency. +function smallestEigenvector(matrix: number[][]): Vec3 { + const work = matrix.map((row) => [...row]); + // Accumulated rotations; its columns are the eigenvectors. + const vectors = [ + [1, 0, 0], + [0, 1, 0], + [0, 0, 1], + ]; + for (let sweep = 0; sweep < 24; sweep++) { + let offDiagonal = 0; + for (let p = 0; p < 3; p++) + for (let q = p + 1; q < 3; q++) offDiagonal += work[p][q] * work[p][q]; + if (offDiagonal < 1e-30) break; + for (let p = 0; p < 3; p++) + for (let q = p + 1; q < 3; q++) { + if (Math.abs(work[p][q]) < 1e-30) continue; + // Standard cyclic-Jacobi rotation (Golub & Van Loan, Matrix + // Computations §8.4): θ = (a_qq − a_pp)/(2·a_pq), and the smaller root + // t = sign(θ)/(|θ| + √(θ²+1)) keeps the rotation angle under 45°. + // θ = 0 is the symmetric case, where either root serves; take +. + const theta = (work[q][q] - work[p][p]) / (2 * work[p][q]); + const tan = + (theta < 0 ? -1 : 1) / + (Math.abs(theta) + Math.sqrt(theta * theta + 1)); + const cos = 1 / Math.sqrt(tan * tan + 1); + const sin = tan * cos; + for (let k = 0; k < 3; k++) { + const kp = work[k][p], + kq = work[k][q]; + work[k][p] = cos * kp - sin * kq; + work[k][q] = sin * kp + cos * kq; + } + for (let k = 0; k < 3; k++) { + const pk = work[p][k], + qk = work[q][k]; + work[p][k] = cos * pk - sin * qk; + work[q][k] = sin * pk + cos * qk; + } + for (let k = 0; k < 3; k++) { + const kp = vectors[k][p], + kq = vectors[k][q]; + vectors[k][p] = cos * kp - sin * kq; + vectors[k][q] = sin * kp + cos * kq; + } + } + } + let best = 0; + for (let i = 1; i < 3; i++) if (work[i][i] < work[best][best]) best = i; + return [vectors[0][best], vectors[1][best], vectors[2][best]]; +} + +// The element boundary faces lying entirely on a picked node set — the same +// membership test the surface loads use, so "the picked surface" means one +// thing across the app. Triangles only: the axis fit needs facet normals, and a +// hex face is split into its two triangles by the caller's winding either way. +function pickedTriangles( + nodeIds: Set, + elements: CouplingElement[], +): [number, number, number][] { + const TET_FACES = [ + [0, 1, 2], + [0, 1, 3], + [0, 2, 3], + [1, 2, 3], + ]; + const HEX_FACES = [ + [0, 1, 2, 3], + [4, 5, 6, 7], + [0, 1, 5, 4], + [1, 2, 6, 5], + [2, 3, 7, 6], + [3, 0, 4, 7], + ]; + const seen = new Set(); + const tris: [number, number, number][] = []; + const add = (a: number, b: number, c: number) => { + const key = [a, b, c].sort((x, y) => x - y).join(","); + if (seen.has(key)) return; + seen.add(key); + tris.push([a, b, c]); + }; + for (const el of elements) { + if (el.type === "CTRIA3") { + if (el.nodeIds.every((n) => nodeIds.has(n))) + add(el.nodeIds[0], el.nodeIds[1], el.nodeIds[2]); + continue; + } + const local = + el.type === "CTETRA" ? TET_FACES : el.type === "CHEXA" ? HEX_FACES : null; + if (!local) continue; + for (const face of local) { + const verts = face.map((i) => el.nodeIds[i]); + if (!verts.every((n) => nodeIds.has(n))) continue; + add(verts[0], verts[1], verts[2]); + if (verts.length === 4) add(verts[0], verts[2], verts[3]); + } + } + return tris; +} + +// Axis of the cylinder a picked surface lies on, or null when it is not one. +// +// Every normal of a cylinder is perpendicular to its axis, so the axis is the +// direction â that minimises Σ A·(n̂·â)² — the smallest eigenvector of the +// area-weighted normal covariance Σ A·n̂n̂ᵀ. Fitting the POSITIONS instead cannot +// do this: for a cylinder of radius R and length L the position covariance has +// eigenvalues L²/12 along the axis and R²/2 across it, so which one is smallest +// flips between a long tube and a short bore — the axis would be found for one +// and lost for the other. +// +// The fit is then CHECKED rather than trusted: the radial distances of the +// picked nodes from the fitted axis must agree to within `tolerance` of their +// mean. A flat face or a fillet passes the eigen-solve just as happily as a bore +// does, and offering "axis end" positions on one would place the reference point +// somewhere with no meaning. +function fitCylinderAxis( + points: Vec3[], + triangles: [number, number, number][], + positionOf: (nodeId: number) => Vec3 | undefined, + tolerance: number, +): { axis: Vec3; centre: Vec3 } | null { + if (points.length < 6 || triangles.length < 2) return null; + + const cov = [ + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + ]; + let totalArea = 0; + for (const [ia, ib, ic] of triangles) { + const cornerA = positionOf(ia), + cornerB = positionOf(ib), + cornerC = positionOf(ic); + if (!cornerA || !cornerB || !cornerC) continue; + const normal = cross(sub(cornerB, cornerA), sub(cornerC, cornerA)); + const twiceArea = Math.hypot(normal[0], normal[1], normal[2]); + if (twiceArea < 1e-20) continue; + const area = 0.5 * twiceArea; + const unit: Vec3 = [ + normal[0] / twiceArea, + normal[1] / twiceArea, + normal[2] / twiceArea, + ]; + totalArea += area; + for (let i = 0; i < 3; i++) + for (let j = 0; j < 3; j++) cov[i][j] += area * unit[i] * unit[j]; + } + if (totalArea < 1e-20) return null; + + const axis = smallestEigenvector(cov); + const len = Math.hypot(axis[0], axis[1], axis[2]); + if (len < 1e-12) return null; + const unitAxis: Vec3 = [axis[0] / len, axis[1] / len, axis[2] / len]; + + const centre = centroidOf(points); + + // Radial spread about the fitted axis — constant on a cylinder, wildly + // uneven on anything else. + const radii = points.map((point) => { + const offset = sub(point, centre); + const along = dot(offset, unitAxis); + return Math.hypot( + offset[0] - along * unitAxis[0], + offset[1] - along * unitAxis[1], + offset[2] - along * unitAxis[2], + ); + }); + const mean = radii.reduce((sum, r) => sum + r, 0) / radii.length; + if (mean < 1e-12) return null; + const variance = + radii.reduce((sum, r) => sum + (r - mean) ** 2, 0) / radii.length; + if (Math.sqrt(variance) / mean > tolerance) return null; + + return { axis: unitAxis, centre }; +} + +// Where a coupling's reference point may sit, given the surface it grips. +// +// Always the centre of the selection — the mean of its nodes, which for a full +// bore is its mid-axis point and for a picked rim is the centre of the ring. +// When the selection is a cylinder (a bolt hole, a bearing seat), the two ends +// of its axis are offered as well: KOF-208's "up and down centre", the positions +// a bolt head and a nut actually occupy. Any other position is typed in as +// coordinates, which is why this list is a starting point and not a constraint. +export function referencePointOptions( + faces: { nodeIds: number[] }[], + nodes: CouplingNode[], + elements: CouplingElement[], + { cylinderTolerance = 0.12 }: { cylinderTolerance?: number } = {}, +): ReferencePointOption[] { + const nodeIds = new Set(); + for (const face of faces) for (const id of face.nodeIds) nodeIds.add(id); + if (nodeIds.size === 0) return []; + + const byId = new Map(nodes.map((n) => [n.id, n])); + const positionOf = (id: number): Vec3 | undefined => { + const n = byId.get(id); + return n ? [n.x, n.y, n.z] : undefined; + }; + const points: Vec3[] = []; + for (const id of nodeIds) { + const position = positionOf(id); + if (position) points.push(position); + } + if (points.length === 0) return []; + + const centroid = centroidOf(points); + const options: ReferencePointOption[] = [ + { label: "Selection centre", point: centroid }, + ]; + + const fit = fitCylinderAxis( + points, + pickedTriangles(nodeIds, elements), + positionOf, + cylinderTolerance, + ); + if (!fit) return options; + + let lo = Infinity, + hi = -Infinity; + for (const point of points) { + const along = dot(sub(point, fit.centre), fit.axis); + if (along < lo) lo = along; + if (along > hi) hi = along; + } + if (!(hi > lo)) return options; + + const at = (along: number): Vec3 => [ + fit.centre[0] + along * fit.axis[0], + fit.centre[1] + along * fit.axis[1], + fit.centre[2] + along * fit.axis[2], + ]; + // Name the two ends by the axis direction rather than "up"/"down": the fitted + // axis has no preferred sign, and a bore drilled along −Z would have its ends + // labelled backwards by any fixed naming. + options.push( + { label: "Axis end (−)", point: at(lo) }, + { label: "Axis end (+)", point: at(hi) }, + ); + return options; +} + +// ── Solver-facing coupling set ──────────────────────────────────────────────── + +// The CSR coupling set solve_coupled takes: coupling k ties `ref[k]` to +// `solid[offsets[k] .. offsets[k+1])`, with `mpc[k]` selecting the kind and +// `dofMask[k]` the DOFs a kinematic coupling ties. +export interface ReferenceCouplingSet { + ref: number[]; + offsets: number[]; + solid: number[]; + mpc: number[]; + dofMask: number[]; +} + +// Build the coupling set from the model's couplings, over pool node indices. +// +// `poolOf` resolves a store node id to its index in the coupled node pool; it +// throws when a node is not in the pool, which is the right outcome — a coupling +// naming a node the solve does not carry is a broken model, not something to +// quietly drop. +export function buildReferenceCouplings( + couplings: CouplingDefinition[], + poolOf: (nodeId: number, context: string) => number, +): ReferenceCouplingSet { + const ref: number[] = []; + const offsets = [0]; + const solid: number[] = []; + const mpc: number[] = []; + const dofMask: number[] = []; + + // A DOF can be eliminated by only one constraint, so a node coupled by two + // KINEMATIC couplings is a modelling error the engine refuses deep inside the + // reduction. Catch it here, where the coupling that caused it can be named. + const kinematicOwner = new Map(); + + for (const coupling of couplings) { + const coupled = coupledNodeIds(coupling); + if (coupled.length === 0) + throw new Error( + `Coupling "${coupling.name}" grips no node — pick the surface it couples to ` + + "its reference point, or delete the coupling.", + ); + if (coupling.kind === "distributing" && coupled.length < 3) + throw new Error( + `Coupling "${coupling.name}" is distributing but grips only ${coupled.length} ` + + "node(s); an RBE3 average needs at least 3. Pick a larger surface, or make " + + "the coupling kinematic.", + ); + const mask = couplingDofMask(coupling.kind, coupling.dofs); + if (mask === 0) + throw new Error( + `Coupling "${coupling.name}" is kinematic but ties no DOF at all — select at ` + + "least one of Ux…Rz.", + ); + + const refPool = poolOf(coupling.refNodeId, `coupling "${coupling.name}"`); + for (const nodeId of coupled) { + if (coupling.kind === "kinematic") { + const owner = kinematicOwner.get(nodeId); + if (owner !== undefined) + throw new Error( + `Node ${nodeId} is gripped by both kinematic couplings "${owner}" and ` + + `"${coupling.name}" — a DOF can be eliminated by only one constraint. ` + + "Make one of them distributing, or pick surfaces that do not overlap.", + ); + kinematicOwner.set(nodeId, coupling.name); + } + solid.push(poolOf(nodeId, `coupling "${coupling.name}" surface`)); + } + ref.push(refPool); + offsets.push(solid.length); + mpc.push(couplingMpcCode(coupling.kind)); + dofMask.push(mask); + } + + return { ref, offsets, solid, mpc, dofMask }; +} + +// Reference-point node ids of a model's couplings — the nodes that carry six +// DOFs without belonging to any element, so the solve path knows to put them in +// the pool and the load path knows a moment on them is a real moment. +export function referencePointIds( + couplings: CouplingDefinition[], +): Set { + return new Set(couplings.map((coupling) => coupling.refNodeId)); +} diff --git a/web/src/lib/facePick.ts b/web/src/lib/facePick.ts index df09409..3747b11 100644 --- a/web/src/lib/facePick.ts +++ b/web/src/lib/facePick.ts @@ -78,7 +78,7 @@ export function sameFaceNodes(a: number[], b: number[]): boolean { export function toggleFaceSelection( current: PickedFace[], picked: PickedFace, - noun: "Face" | "Edge" = "Face", + noun: SelectionNoun = "Face", ): PickedFace[] { const existingIdx = current.findIndex((f) => sameFaceNodes(f.nodeIds, picked.nodeIds), @@ -89,10 +89,34 @@ export function toggleFaceSelection( : [...current, picked]; return next.map((f, i) => ({ ...f, - label: `${noun} ${i + 1} (${f.nodeIds.length} nodes)`, + label: selectionLabel(noun, i + 1, f.nodeIds.length), })); } +// What a click selects. "face" is a surface region, "edge" the boundary +// polyline near the click, "point" the single nearest node of the clicked +// facet. +export type PickGeometry = "face" | "edge" | "point"; + +export type SelectionNoun = "Face" | "Edge" | "Node"; + +export const SELECTION_NOUN: Record = { + face: "Face", + edge: "Edge", + point: "Node", +}; + +// "Face 1 (37 nodes)" / "Node 2 (1 node)". Written in one place because both +// the viewport's live labels and the panel's committed entries use it, and a +// selection that reads differently in the two is the same selection twice. +export function selectionLabel( + noun: SelectionNoun, + index: number, + nNodes: number, +): string { + return `${noun} ${index} (${nNodes} node${nNodes === 1 ? "" : "s"})`; +} + export function buildEdgeToTris(triangles: Tri[]): Map { const edgeToTris = new Map(); for (let i = 0; i < triangles.length; i++) { @@ -437,3 +461,67 @@ export function pickEdgeNodeIds( } return nodeIds; } + +// ── Point picking ───────────────────────────────────────────────────────────── + +/** + * The reference point nearest a click, and how far away it was. + * + * A reference point is not part of any surface, so the ray never reports it — + * the click always lands on the mesh behind it. Whether the user meant the point + * or the node behind it is therefore a DISTANCE question, decided here rather + * than left to the depth ordering of two overlapping meshes: a marker drawn on + * the face it couples sits fractionally behind that face, so the surface always + * wins the raycast even when the point is what was clicked. + */ +export function nearestReferencePoint( + clickPoint: Vec3, + points: { nodeId: number; point: Vec3 }[], +): { nodeId: number; distSq: number } | null { + let nearest: { nodeId: number; distSq: number } | null = null; + for (const { nodeId, point } of points) { + const distSq = + (clickPoint[0] - point[0]) ** 2 + + (clickPoint[1] - point[1]) ** 2 + + (clickPoint[2] - point[2]) ** 2; + if (!nearest || distSq < nearest.distSq) nearest = { nodeId, distSq }; + } + return nearest; +} + +/** + * Pick the single node nearest the click, among the vertices of the clicked + * facet. + * + * Restricting the candidates to the clicked triangle is what makes this exact: + * a global nearest-node search over the whole mesh would happily return a node + * on the far side of the model that happens to project closer in space, and the + * user would have selected something they cannot see. The ray already told us + * which facet was hit; the node they meant is one of its three corners. + * + * Returns a one-element set, the same node-id-set shape the face and edge picks + * return, so group creation, storage and the viewport need no special case. + */ +export function pickPointNodeId( + clickPoint: Vec3, + startTriIdx: number, + topo: BoundaryMeshTopo, + getPos: (id: number) => Vec3, +): Set { + const triangle = topo.triangles[startTriIdx]; + if (!triangle) return new Set(); + let nearest = triangle[0]; + let best = Infinity; + for (const nodeId of triangle) { + const pos = getPos(nodeId); + const distSq = + (clickPoint[0] - pos[0]) ** 2 + + (clickPoint[1] - pos[1]) ** 2 + + (clickPoint[2] - pos[2]) ** 2; + if (distSq < best) { + best = distSq; + nearest = nodeId; + } + } + return new Set([nearest]); +} diff --git a/web/src/lib/shellize.ts b/web/src/lib/shellize.ts index ea30d8e..613cc8f 100644 --- a/web/src/lib/shellize.ts +++ b/web/src/lib/shellize.ts @@ -667,16 +667,36 @@ export function extractThinWallShells( // CSR-style coupling set. `mpc[k]` selects the coupling kind for reference k: // 1 ⇒ relaxed shell-to-solid MPC (rigid translation tie + relaxed rotation, for a -// continuous-material seam), 0/absent ⇒ distributing RBE3 (for a gapped interface). +// continuous-material seam), 2 ⇒ kinematic RBE2 (the coupled nodes follow the +// reference point rigidly), 0/absent ⇒ distributing RBE3 (for a gapped +// interface). `dofMask[k]` selects which of a kinematic coupling's six DOFs are +// tied (bits 0..5); absent ⇒ all six, which is what every coupling this module +// derives itself wants. export interface CouplingSet { ref: number[]; offsets: number[]; solid: number[]; mpc?: number[]; // per ref; length matches ref when present + dofMask?: number[]; // per ref; length matches ref when present +} + +// All six DOFs of a coupling — the engine's kAllDofs. +const ALL_DOF_MASK = 0x3f; + +// Per-reference coupling kinds / DOF masks of a set, defaulted for the sets +// built here (which are all distributing or relaxed-MPC, and tie all six DOFs). +export function couplingMpcCodes(set: CouplingSet): number[] { + if (set.mpc) return set.mpc; + return set.ref.map(() => 0); +} + +export function couplingDofMasks(set: CouplingSet): number[] { + if (set.dofMask) return set.dofMask; + return set.ref.map(() => ALL_DOF_MASK); } export interface CoupledModel { - pool: number[]; // 3·nPool (solid nodes then shell nodes) + pool: number[]; // 3·nPool (solid nodes, then shell nodes, then reference points) tets: number[]; // 4·nSolidTets over pool tetBody: number[]; // body id per solid tet (for per-body PSOLID labelling) triangles: number[]; // 3·nShellTris over pool @@ -684,6 +704,7 @@ export interface CoupledModel { coupling: CouplingSet; solidPool: Map; // original vertex index → pool index (solid) shellPool: number[]; // local shell index → pool index + refPool: Map; // reference-point vertex index → pool index } // Boundary of a tet mesh: the faces used by exactly one element, as a flat @@ -889,19 +910,22 @@ function autoDetectCouplings( return { ref, offsets, solid }; } -// Concatenate two CSR coupling sets (shell↔solid and solid↔solid) into one, -// preserving each set's per-reference MPC flags (missing ⇒ distributing). -function concatCouplings(a: CouplingSet, b: CouplingSet): CouplingSet { +// Concatenate two CSR coupling sets (shell↔solid, solid↔solid, and the +// reference-point couplings the user declared) into one, preserving each set's +// per-reference kind and DOF mask. +export function concatCouplings(a: CouplingSet, b: CouplingSet): CouplingSet { const ref = [...a.ref, ...b.ref]; const solid = [...a.solid, ...b.solid]; const offsets = [...a.offsets]; const base = a.solid.length; for (let k = 1; k < b.offsets.length; k++) offsets.push(base + b.offsets[k]); - const mpc = [ - ...(a.mpc ?? a.ref.map(() => 0)), - ...(b.mpc ?? b.ref.map(() => 0)), - ]; - return { ref, offsets, solid, mpc }; + return { + ref, + offsets, + solid, + mpc: [...couplingMpcCodes(a), ...couplingMpcCodes(b)], + dofMask: [...couplingDofMasks(a), ...couplingDofMasks(b)], + }; } // One tie connection, as the coupled builders need it: the two surfaces the user @@ -1055,20 +1079,27 @@ export function dropCouplingsOnFixedNodes( ): CouplingSet { const fixedNodes = new Set(); for (const d of fixedDofs) fixedNodes.add(Math.floor(d / 6)); + const kinds = couplingMpcCodes(coupling); + const masks = couplingDofMasks(coupling); const ref: number[] = [], offsets = [0], solid: number[] = [], - mpc: number[] = []; + mpc: number[] = [], + dofMask: number[] = []; for (let k = 0; k < coupling.ref.length; k++) { - if (fixedNodes.has(coupling.ref[k])) continue; + // A KINEMATIC coupling's reference point is INDEPENDENT — fixing it is the + // whole point of a bolted or clamped reference point, not a conflict — so + // only the couplings whose reference is dependent (distributing, relaxed + // MPC) are dropped here. + if (kinds[k] !== 2 && fixedNodes.has(coupling.ref[k])) continue; ref.push(coupling.ref[k]); for (let i = coupling.offsets[k]; i < coupling.offsets[k + 1]; i++) solid.push(coupling.solid[i]); offsets.push(solid.length); - // eslint-disable-next-line kofem/no-silent-fallback -- mpc is an optional coupling-kind flag; absent means the distributing (RBE3) default - mpc.push(coupling.mpc?.[k] ?? 0); + mpc.push(kinds[k]); + dofMask.push(masks[k]); } - return { ref, offsets, solid, mpc }; + return { ref, offsets, solid, mpc, dofMask }; } // Assemble the coupled node pool (solid nodes then shell nodes), remap tets and @@ -1084,6 +1115,7 @@ export function buildCoupledModel( couplingRadius, maxCoupledNodes = 16, ties = [], + referencePoints = [], }: { // How far off the retained solid's boundary a shell node may sit and still // count as seam. Defaults to seamTolerance() of the model's own scales. @@ -1096,6 +1128,9 @@ export function buildCoupledModel( // these — nothing is inferred from the geometry — so an assembly with no // connection keeps its bodies apart, which is what the empty default means. ties?: TieSurfaces[]; + // Vertex indices of the surface-to-point couplings' reference points. They + // belong to no tet, so nothing else puts them in the pool. + referencePoints?: number[]; } = {}, ): CoupledModel { // Solid tets = the other bodies plus the shelled body's non-wall (base) tets; @@ -1132,6 +1167,11 @@ export function buildCoupledModel( shells.shellVerts[3 * s + 2], ), ); + // Reference points last, so isShellPoolIndex can keep telling the three + // groups apart by index range. + const refPool = new Map(); + for (const vi of referencePoints) + if (!refPool.has(vi)) refPool.set(vi, addPool(...pt(m.V, vi))); const tets = solidTets.map((n) => { const poolIndex = solidPool.get(n); @@ -1184,6 +1224,7 @@ export function buildCoupledModel( thicknesses: shells.shellThk, solidPool, shellPool, + refPool, // The shell↔solid seam is continuous material (a thin wall idealised as shell, // tied back to its retained solid) — not a tie connection between parts, and // never something the user declares — so it uses the relaxed MPC coupling @@ -1226,6 +1267,7 @@ export function buildExplicitCoupledModel( couplingRadius, maxCoupledNodes = 16, ties = [], + referencePoints = [], }: { seamTolerance?: number; couplingRadius?: number; @@ -1233,6 +1275,10 @@ export function buildExplicitCoupledModel( // The model's tie connections — the only thing that joins distinct solid // bodies here (see buildCoupledModel). ties?: TieSurfaces[]; + // Store vertex indices of the surface-to-point couplings' REFERENCE POINTS. + // They belong to no element, so nothing else would put them in the pool — + // and without a pool node the coupling has no reference to tie to. + referencePoints?: number[]; } = {}, ): ExplicitCoupledModel { const pool: number[] = []; @@ -1262,6 +1308,10 @@ export function buildExplicitCoupledModel( shellPoolIndex.add(pi); return pi; }); + // Reference points last. They carry no element stiffness; the engine gives a + // coupling reference its six DOFs and leaves its rotations free (shell_core: + // is_coupling_ref), so they need nothing here beyond a place in the pool. + for (const vi of referencePoints) addVertex(vi); const ppt = (i: number): [number, number, number] => [ pool[3 * i], @@ -1323,13 +1373,18 @@ export function buildExplicitCoupledModel( }; } -// A shell node is "solid" in the pool iff its index is < solidPool.size; shell -// nodes are appended after the solid nodes. +// The pool is built in three blocks — solid nodes, then shell mid-surface +// nodes, then reference points — so a pool index says which it is. Only the +// middle block carries shell (6-DOF) element stiffness; a reference point has +// six DOFs too, but they come from its coupling, not from a facet. export function isShellPoolIndex( model: CoupledModel, poolIndex: number, ): boolean { - return poolIndex >= model.solidPool.size; + return ( + poolIndex >= model.solidPool.size && + poolIndex < model.solidPool.size + model.shellPool.length + ); } // Nearest shell-mid-surface pool node to a point — used to map the shelled diff --git a/web/src/store/boundarySlice.ts b/web/src/store/boundarySlice.ts index 7475897..8a9e352 100644 --- a/web/src/store/boundarySlice.ts +++ b/web/src/store/boundarySlice.ts @@ -8,6 +8,9 @@ import type { SliceCreator } from "./modelStore"; import type { Element, Node } from "./geometrySlice"; import type { TieDefinition, TieExtent } from "../lib/tie"; +import type { PickGeometry } from "../lib/facePick"; +import type { CouplingDefinition, CouplingKind } from "../lib/coupling"; +import { ALL_DOFS, referencePointIds } from "../lib/coupling"; import { momentToNodalForces } from "./momentLoad"; export interface Constraint { @@ -35,6 +38,12 @@ export interface BcFaceEntry { id: number; label: string; // e.g. "Face 1" nodeIds: number[]; + // What was picked. Only "point" changes how the solver applies the group: a + // single node spans no element face, so a load on it cannot be integrated as + // a traction and is applied at the node instead (rebuildLoads). Absent on + // entries saved before point picking existed, which were faces or edges — + // both integrable — so absence reads correctly as "not a point". + geometry?: PickGeometry; } export interface NamedBcGroup { @@ -81,9 +90,32 @@ export interface TieGroup extends TieDefinition { searchDistance: number; } +// A surface-to-point coupling — the picked surface, the reference point it is +// idealised to, and how the two are tied (see lib/coupling.ts). +// +// The reference point is a real node of the model: it is appended to `nodes` +// when the coupling is created and removed with it. That is what a reference +// point IS in a solver deck, and it means a BC or a load reaches it through +// exactly the machinery every other node uses — `refNodeId` in a face entry, +// nothing special anywhere downstream. The invariant it rests on: a node that +// belongs to no element exists only because a coupling created it, so deleting +// the coupling must delete the node (deleteCouplingGroup does). +export interface CouplingGroup extends CouplingDefinition { + id: number; + name: string; // e.g. "Coupling1" + faces: BcFaceEntry[]; + // The reference point's position, mirrored from the node so the panel can + // edit it without reaching into the mesh. + point: [number, number, number]; +} + // Which surface of a tie a pick session is currently filling. export type TieSide = "a" | "b"; +// What a pick session is being used to define. A coupling picks one surface, +// like a BC or a load; a tie picks two (pickTieSide selects which). +export type PickMode = "bc" | "load" | "tie" | "coupling"; + // Default search distance (mm) offered for a region tie. A contact patch is a // property of the assembly, so there is no right number — this is only the // starting value in the form, which the user edits before applying. @@ -242,21 +274,91 @@ function loadedFaces( return edges; } +// Mint a committed face entry from a picked one, taking the next id. Written +// once because every group type creates entries the same way, and a copy that +// forgets to carry a field — `geometry` is the one that matters, since it is +// what routes a load to its nodes rather than to a surface integral — fails +// only in the solve, far from here. +function faceEntry( + face: Omit, + nextId: () => number, +): BcFaceEntry { + return { ...face, id: nextId() }; +} + +// Whether a selection is applied AT ITS NODES rather than integrated over a +// surface. Two ways to be one, and neither spans an element face: +// +// a coupling's reference point — a node belonging to no element at all +// a point pick — a single mesh node, which spans no face +// +// The single definition rebuildLoads and rebuildSurfaceLoads both consult to +// decide which of them owns a face. Written once on purpose: a load that both +// claim is applied twice, and a load neither claims disappears without a word. +function isNodalFace( + face: { nodeIds: number[]; geometry?: PickGeometry }, + refPoints: Set, +): boolean { + if (face.geometry === "point") return true; + return ( + face.nodeIds.length > 0 && face.nodeIds.every((id) => refPoints.has(id)) + ); +} + export function rebuildLoads( loadGroups: NamedLoadGroup[], nodes: Node[], + couplingGroups: CouplingGroup[] = [], ): Load[] { const nodeById = new Map(); for (const n of nodes) nodeById.set(n.id, n); + const refPoints = referencePointIds(couplingGroups); const result: Load[] = []; for (const g of loadGroups) { - // Force and pressure loads are applied as work-equivalent surface tractions - // (rebuildSurfaceLoads), not lumped nodal forces — they are skipped here. - if (loadKind(g) !== "moment") continue; - result.push( - ...momentToNodalForces(loadComponents(g), g.faces, nodeById, g.name), + const kind = loadKind(g); + if (kind === "pressure") continue; // no vector, and a point has no area + const vector = loadComponents(g); + + // A load on a single NODE — a picked point, or a coupling's reference point + // — is applied there directly, whichever kind it is, because neither of the + // usual routes can carry it: + // + // force — rebuildSurfaceLoads integrates a traction over the boundary + // faces lying on the selection, and a lone node spans none, so + // the load would be integrated over nothing and vanish. + // moment — momentToNodalForces turns a couple into a ring of tangential + // forces about the face centroid, and a single node IS its own + // centroid, so every lever arm is zero and the couple vanishes. + // + // A reference point carries six real DOFs (it is a coupling reference — see + // shell_core's is_coupling_ref), so it takes both the force and the couple, + // and its coupling spreads them over the surface it grips. An ordinary mesh + // node has only the three translations; a moment applied there reaches the + // engine as a rotational DOF the solid solve does not carry, so it is + // rejected at the panel rather than accepted and dropped here. + // + // The group's vector is its TOTAL, so several nodes in one group share it, + // exactly as a surface selection shares a traction. + const pointFaces = new Set( + g.faces.filter((face) => isNodalFace(face, refPoints)), ); + const pointNodeIds = [...pointFaces].flatMap((face) => face.nodeIds); + for (const nodeId of pointNodeIds) + for (let axis = 0; axis < 3; axis++) + if (vector[axis] !== 0) + result.push({ + nodeId, + dof: (kind === "moment" ? 3 : 0) + axis, + value: vector[axis] / pointNodeIds.length, + }); + + // Faces of ordinary mesh nodes keep their existing routes: a force is a + // work-equivalent surface traction (rebuildSurfaceLoads, not here), a moment + // becomes equivalent nodal forces. + if (kind !== "moment") continue; + const meshFaces = g.faces.filter((face) => !pointFaces.has(face)); + result.push(...momentToNodalForces(vector, meshFaces, nodeById, g.name)); } return result; } @@ -272,12 +374,21 @@ export function rebuildLoads( export function rebuildSurfaceLoads( loadGroups: NamedLoadGroup[], elements: Element[], + couplingGroups: CouplingGroup[] = [], ): SurfaceLoad[] { + const refPoints = referencePointIds(couplingGroups); const result: SurfaceLoad[] = []; for (const g of loadGroups) { const kind = loadKind(g); if (kind === "moment") continue; // moments stay as equivalent point loads for (const f of g.faces) { + // Nodal selections are rebuildLoads' half of the group: a reference point + // belongs to no element and a picked point spans no face, so there is no + // surface to integrate over. The two builders partition a group's faces by + // the SAME test on purpose — leaving it to loadedFaces returning nothing + // would make the split an accident of that function, and a change there + // would either double-apply the load or drop it without a word. + if (isNodalFace(f, refPoints)) continue; const faces = loadedFaces(f, elements); if (faces.length === 0) continue; if (kind === "pressure") { @@ -299,22 +410,25 @@ export interface BoundarySlice { // work-equivalent surface tractions (force/pressure groups) handed to the solver surfaceLoads: SurfaceLoad[]; - // Named BC / Load / Tie groups (primary source of truth for constraints, - // loads and bonded connections) + // Named BC / Load / Tie / Coupling groups (primary source of truth for + // constraints, loads, bonded connections and surface-to-point couplings) bcGroups: NamedBcGroup[]; loadGroups: NamedLoadGroup[]; tieGroups: TieGroup[]; + couplingGroups: CouplingGroup[]; nextBcGroupId: number; nextLoadGroupId: number; nextTieGroupId: number; + nextCouplingGroupId: number; nextFaceEntryId: number; - pickMode: "bc" | "load" | "tie" | null; + pickMode: PickMode | null; pickTargetGroupId: number | null; // null = creating new group; id = adding to existing - // Whether a click selects a surface region ("face") or a boundary polyline - // ("edge"). Edge picking is the only way to grab the rim of a flat shell, - // whose whole sheet is a single face-pick region. Reset to "face" on exit. - pickGeometry: "face" | "edge"; + // What a click selects: a surface region ("face"), the boundary polyline near + // the click ("edge" — the only way to grab the rim of a flat shell, whose + // whole sheet is a single face-pick region), or the single nearest node + // ("point"). Reset to "face" on exit. + pickGeometry: PickGeometry; selectedFace: FaceSelection | null; pendingFaces: FaceSelection[]; // faces accumulated via shift-click within a pick session @@ -325,15 +439,24 @@ export interface BoundarySlice { pickTieSide: TieSide; tieDraft: { a: FaceSelection[]; b: FaceSelection[] }; + // The reference point currently being placed, and the nodes it would grip — + // the coupling as it stands in a form that has not been applied yet. A + // reference point is a position in space with nothing in the mesh to anchor it + // to, so typing coordinates without seeing them is guesswork; the viewport + // draws this the way it draws a committed coupling. Null whenever no coupling + // form is open. The form that owns it supplies both fields, rather than the + // viewport inferring the nodes from whichever session it thinks is live. + couplingDraft: { point: [number, number, number]; nodeIds: number[] } | null; + // Pick mode / face selection - setPickMode( - mode: "bc" | "load" | "tie" | null, - targetGroupId?: number | null, - ): void; - setPickGeometry(geometry: "face" | "edge"): void; + setPickMode(mode: PickMode | null, targetGroupId?: number | null): void; + setPickGeometry(geometry: PickGeometry): void; setSelectedFace(face: FaceSelection | null): void; setPendingFaces(faces: FaceSelection[]): void; setPickTieSide(side: TieSide): void; + setCouplingDraft( + draft: { point: [number, number, number]; nodeIds: number[] } | null, + ): void; // BC group actions createBcGroup( @@ -382,6 +505,25 @@ export interface BoundarySlice { removeFaceFromTieGroup(groupId: number, side: TieSide, faceId: number): void; deleteTieGroup(id: number): void; clearTies(): void; + + // Surface-to-point coupling actions. Creating one also creates its reference + // point node; deleting one removes that node again. + createCouplingGroup( + faces: Omit[], + point: [number, number, number], + kind: CouplingKind, + dofs: number[], + ): void; + addFaceToCouplingGroup(groupId: number, face: Omit): void; + updateCouplingGroup( + id: number, + kind: CouplingKind, + dofs: number[], + point: [number, number, number], + ): void; + removeFaceFromCouplingGroup(groupId: number, faceId: number): void; + deleteCouplingGroup(id: number): void; + clearCouplings(): void; } export const createBoundarySlice: SliceCreator = (set) => ({ @@ -391,9 +533,11 @@ export const createBoundarySlice: SliceCreator = (set) => ({ bcGroups: [], loadGroups: [], tieGroups: [], + couplingGroups: [], nextBcGroupId: 1, nextLoadGroupId: 1, nextTieGroupId: 1, + nextCouplingGroupId: 1, nextFaceEntryId: 1, pickMode: null, pickTargetGroupId: null, @@ -402,17 +546,18 @@ export const createBoundarySlice: SliceCreator = (set) => ({ pendingFaces: [], pickTieSide: "a", tieDraft: { a: [], b: [] }, + couplingDraft: null, // Pick mode / face selection - setPickMode: ( - mode: "bc" | "load" | "tie" | null, - targetGroupId: number | null = null, - ) => + setPickMode: (mode: PickMode | null, targetGroupId: number | null = null) => set((s) => { s.pickMode = mode; s.pickTargetGroupId = mode !== null ? (targetGroupId ?? null) : null; s.pickTieSide = "a"; s.tieDraft = { a: [], b: [] }; + // Leaving the pick session abandons the coupling being placed with it, so + // its preview must not outlive the form that owned it. + s.couplingDraft = null; if (mode === null) { s.selectedFace = null; s.pendingFaces = []; @@ -420,6 +565,11 @@ export const createBoundarySlice: SliceCreator = (set) => ({ } }), + setCouplingDraft: (draft) => + set((s) => { + s.couplingDraft = draft; + }), + // Park the side being picked and bring the other one into the live session, // so a tie's two surfaces are filled by the same face-pick machinery. setPickTieSide: (side: TieSide) => @@ -436,7 +586,7 @@ export const createBoundarySlice: SliceCreator = (set) => ({ s.pickTieSide = side; }), - setPickGeometry: (geometry: "face" | "edge") => + setPickGeometry: (geometry: PickGeometry) => set((s) => { s.pickGeometry = geometry; }), @@ -458,11 +608,9 @@ export const createBoundarySlice: SliceCreator = (set) => ({ value: number, ) => set((s) => { - const faceEntries = faces.map((f) => ({ - id: s.nextFaceEntryId++, - label: f.label, - nodeIds: f.nodeIds, - })); + const faceEntries = faces.map((f) => + faceEntry(f, () => s.nextFaceEntryId++), + ); s.bcGroups.push({ id: s.nextBcGroupId, name: `BC${s.nextBcGroupId}`, @@ -479,12 +627,7 @@ export const createBoundarySlice: SliceCreator = (set) => ({ set((s) => { const group = s.bcGroups.find((g) => g.id === groupId); if (!group) return; - const faceId = s.nextFaceEntryId++; - group.faces.push({ - id: faceId, - label: face.label, - nodeIds: face.nodeIds, - }); + group.faces.push(faceEntry(face, () => s.nextFaceEntryId++)); s.constraints = rebuildConstraints(s.bcGroups); s.result = null; }), @@ -533,11 +676,9 @@ export const createBoundarySlice: SliceCreator = (set) => ({ components?: [number, number, number], ) => set((s) => { - const faceEntries = faces.map((f) => ({ - id: s.nextFaceEntryId++, - label: f.label, - nodeIds: f.nodeIds, - })); + const faceEntries = faces.map((f) => + faceEntry(f, () => s.nextFaceEntryId++), + ); // A componentwise force/moment carries its vector in `components` (the // source of truth); dof/totalForce are derived as a legacy summary. const { dof, totalForce } = components @@ -553,8 +694,12 @@ export const createBoundarySlice: SliceCreator = (set) => ({ kind, }); s.nextLoadGroupId++; - s.loads = rebuildLoads(s.loadGroups, s.nodes); - s.surfaceLoads = rebuildSurfaceLoads(s.loadGroups, s.elements); + s.loads = rebuildLoads(s.loadGroups, s.nodes, s.couplingGroups); + s.surfaceLoads = rebuildSurfaceLoads( + s.loadGroups, + s.elements, + s.couplingGroups, + ); s.result = null; }), @@ -562,14 +707,13 @@ export const createBoundarySlice: SliceCreator = (set) => ({ set((s) => { const group = s.loadGroups.find((g) => g.id === groupId); if (!group) return; - const faceId = s.nextFaceEntryId++; - group.faces.push({ - id: faceId, - label: face.label, - nodeIds: face.nodeIds, - }); - s.loads = rebuildLoads(s.loadGroups, s.nodes); - s.surfaceLoads = rebuildSurfaceLoads(s.loadGroups, s.elements); + group.faces.push(faceEntry(face, () => s.nextFaceEntryId++)); + s.loads = rebuildLoads(s.loadGroups, s.nodes, s.couplingGroups); + s.surfaceLoads = rebuildSurfaceLoads( + s.loadGroups, + s.elements, + s.couplingGroups, + ); s.result = null; }), @@ -594,8 +738,12 @@ export const createBoundarySlice: SliceCreator = (set) => ({ group.totalForce = totalForce; if (components) group.components = components; else delete group.components; - s.loads = rebuildLoads(s.loadGroups, s.nodes); - s.surfaceLoads = rebuildSurfaceLoads(s.loadGroups, s.elements); + s.loads = rebuildLoads(s.loadGroups, s.nodes, s.couplingGroups); + s.surfaceLoads = rebuildSurfaceLoads( + s.loadGroups, + s.elements, + s.couplingGroups, + ); s.result = null; }), @@ -606,16 +754,24 @@ export const createBoundarySlice: SliceCreator = (set) => ({ group.faces = group.faces.filter((f) => f.id !== faceId); if (group.faces.length === 0) s.loadGroups = s.loadGroups.filter((g) => g.id !== groupId); - s.loads = rebuildLoads(s.loadGroups, s.nodes); - s.surfaceLoads = rebuildSurfaceLoads(s.loadGroups, s.elements); + s.loads = rebuildLoads(s.loadGroups, s.nodes, s.couplingGroups); + s.surfaceLoads = rebuildSurfaceLoads( + s.loadGroups, + s.elements, + s.couplingGroups, + ); s.result = null; }), deleteLoadGroup: (id: number) => set((s) => { s.loadGroups = s.loadGroups.filter((g) => g.id !== id); - s.loads = rebuildLoads(s.loadGroups, s.nodes); - s.surfaceLoads = rebuildSurfaceLoads(s.loadGroups, s.elements); + s.loads = rebuildLoads(s.loadGroups, s.nodes, s.couplingGroups); + s.surfaceLoads = rebuildSurfaceLoads( + s.loadGroups, + s.elements, + s.couplingGroups, + ); s.result = null; }), @@ -638,11 +794,7 @@ export const createBoundarySlice: SliceCreator = (set) => ({ ) => set((s) => { const entries = (faces: Omit[]) => - faces.map((face) => ({ - id: s.nextFaceEntryId++, - label: face.label, - nodeIds: face.nodeIds, - })); + faces.map((face) => faceEntry(face, () => s.nextFaceEntryId++)); s.tieGroups.push({ id: s.nextTieGroupId, name: `Tie${s.nextTieGroupId}`, @@ -663,11 +815,7 @@ export const createBoundarySlice: SliceCreator = (set) => ({ set((s) => { const group = s.tieGroups.find((tie) => tie.id === groupId); if (!group) return; - tieFaces(group, side).push({ - id: s.nextFaceEntryId++, - label: face.label, - nodeIds: face.nodeIds, - }); + tieFaces(group, side).push(faceEntry(face, () => s.nextFaceEntryId++)); s.result = null; }), @@ -706,4 +854,129 @@ export const createBoundarySlice: SliceCreator = (set) => ({ s.tieGroups = []; s.result = null; }), + + // Surface-to-point coupling actions. The reference point is a NODE, created + // and destroyed with its coupling — see CouplingGroup. + createCouplingGroup: ( + faces: Omit[], + point: [number, number, number], + kind: CouplingKind, + dofs: number[], + ) => + set((s) => { + const refNodeId = + s.nodes.reduce((highest, node) => Math.max(highest, node.id), -1) + 1; + s.nodes.push({ id: refNodeId, x: point[0], y: point[1], z: point[2] }); + s.couplingGroups.push({ + id: s.nextCouplingGroupId, + name: `Coupling${s.nextCouplingGroupId}`, + kind, + dofs: kind === "kinematic" ? dofs : ALL_DOFS, + refNodeId, + point, + faces: faces.map((face) => faceEntry(face, () => s.nextFaceEntryId++)), + }); + s.nextCouplingGroupId++; + s.result = null; + }), + + addFaceToCouplingGroup: (groupId: number, face: Omit) => + set((s) => { + const group = s.couplingGroups.find((c) => c.id === groupId); + if (!group) return; + group.faces.push(faceEntry(face, () => s.nextFaceEntryId++)); + s.result = null; + }), + + updateCouplingGroup: ( + id: number, + kind: CouplingKind, + dofs: number[], + point: [number, number, number], + ) => + set((s) => { + const group = s.couplingGroups.find((c) => c.id === id); + if (!group) return; + group.kind = kind; + // A distributing coupling ties all six DOFs of its reference point by + // construction — the mask is a kinematic-only control, so storing a + // partial one would describe a constraint the solver does not apply. + group.dofs = kind === "kinematic" ? dofs : ALL_DOFS; + group.point = point; + const refNode = s.nodes.find((node) => node.id === group.refNodeId); + if (refNode) { + refNode.x = point[0]; + refNode.y = point[1]; + refNode.z = point[2]; + } + // Moving the point changes the lever arms of a moment applied to it. + s.loads = rebuildLoads(s.loadGroups, s.nodes, s.couplingGroups); + s.result = null; + }), + + // A coupling with no surface grips nothing, so losing its last face removes + // it — the same rule a BC group follows. + removeFaceFromCouplingGroup: (groupId: number, faceId: number) => + set((s) => { + const group = s.couplingGroups.find((c) => c.id === groupId); + if (!group) return; + group.faces = group.faces.filter((face) => face.id !== faceId); + if (group.faces.length === 0) removeCoupling(s, groupId); + s.result = null; + }), + + deleteCouplingGroup: (id: number) => + set((s) => { + removeCoupling(s, id); + s.result = null; + }), + + clearCouplings: () => + set((s) => { + for (const group of [...s.couplingGroups]) removeCoupling(s, group.id); + s.result = null; + }), }); + +// Drop a coupling and the reference point node it owns, together with the BCs +// and loads that were applied to that point. Leaving them behind would leave +// constraints on a node that no longer exists — and, because a reference point +// belongs to no element, a free node in the mesh that the all-solid solve would +// assemble into a singular system. +function removeCoupling( + s: { + couplingGroups: CouplingGroup[]; + bcGroups: NamedBcGroup[]; + loadGroups: NamedLoadGroup[]; + nodes: Node[]; + elements: Element[]; + constraints: Constraint[]; + loads: Load[]; + surfaceLoads: SurfaceLoad[]; + }, + id: number, +): void { + const group = s.couplingGroups.find((c) => c.id === id); + if (!group) return; + s.couplingGroups = s.couplingGroups.filter((c) => c.id !== id); + s.nodes = s.nodes.filter((node) => node.id !== group.refNodeId); + + const withoutPoint = (groups: T[]): T[] => + groups + .map((g) => ({ + ...g, + faces: g.faces.filter( + (face) => !face.nodeIds.includes(group.refNodeId), + ), + })) + .filter((g) => g.faces.length > 0); + s.bcGroups = withoutPoint(s.bcGroups); + s.loadGroups = withoutPoint(s.loadGroups); + s.constraints = rebuildConstraints(s.bcGroups); + s.loads = rebuildLoads(s.loadGroups, s.nodes, s.couplingGroups); + s.surfaceLoads = rebuildSurfaceLoads( + s.loadGroups, + s.elements, + s.couplingGroups, + ); +} diff --git a/web/src/store/geometrySlice.ts b/web/src/store/geometrySlice.ts index abbba21..d9370ce 100644 --- a/web/src/store/geometrySlice.ts +++ b/web/src/store/geometrySlice.ts @@ -288,12 +288,14 @@ export const createGeometrySlice: SliceCreator = (set) => ({ s.bcGroups = []; s.loadGroups = []; s.tieGroups = []; + s.couplingGroups = []; s.constraints = []; s.loads = []; s.surfaceLoads = []; s.nextBcGroupId = 1; s.nextLoadGroupId = 1; s.nextTieGroupId = 1; + s.nextCouplingGroupId = 1; s.nextFaceEntryId = 1; s.result = null; if (tessellation) { @@ -337,12 +339,14 @@ export const createGeometrySlice: SliceCreator = (set) => ({ s.bcGroups = []; s.loadGroups = []; s.tieGroups = []; + s.couplingGroups = []; s.constraints = []; s.loads = []; s.surfaceLoads = []; s.nextBcGroupId = 1; s.nextLoadGroupId = 1; s.nextTieGroupId = 1; + s.nextCouplingGroupId = 1; s.nextFaceEntryId = 1; s.result = null; s.selectedFace = null; @@ -351,6 +355,7 @@ export const createGeometrySlice: SliceCreator = (set) => ({ s.pickTargetGroupId = null; s.pickTieSide = "a"; s.tieDraft = { a: [], b: [] }; + s.couplingDraft = null; s.modelName = name; s.viewRepr = "surface"; s.fitViewTrigger++; diff --git a/web/src/store/modelStore.ts b/web/src/store/modelStore.ts index 0e86669..10ac007 100644 --- a/web/src/store/modelStore.ts +++ b/web/src/store/modelStore.ts @@ -58,6 +58,8 @@ export type { SurfaceLoad, TieGroup, TieSide, + PickMode, + CouplingGroup, } from "./boundarySlice"; export { loadKind, @@ -66,6 +68,9 @@ export { DEFAULT_TIE_DISTANCE, } from "./boundarySlice"; export type { TieExtent } from "../lib/tie"; +export type { PickGeometry } from "../lib/facePick"; +export type { CouplingKind, ReferencePointOption } from "../lib/coupling"; +export { ALL_DOFS, referencePointOptions } from "../lib/coupling"; export type { SolverResult, ResultType, @@ -117,12 +122,18 @@ const createAnalysisActions: SliceCreator = (set) => ({ s.bcGroups = a.bcGroups; s.loadGroups = a.loadGroups; s.tieGroups = a.tieGroups; + s.couplingGroups = a.couplingGroups; s.constraints = rebuildConstraints(a.bcGroups); - s.loads = rebuildLoads(a.loadGroups, s.nodes); - s.surfaceLoads = rebuildSurfaceLoads(a.loadGroups, a.elements); + s.loads = rebuildLoads(a.loadGroups, s.nodes, a.couplingGroups); + s.surfaceLoads = rebuildSurfaceLoads( + a.loadGroups, + a.elements, + a.couplingGroups, + ); s.nextBcGroupId = a.nextBcGroupId; s.nextLoadGroupId = a.nextLoadGroupId; s.nextTieGroupId = a.nextTieGroupId; + s.nextCouplingGroupId = a.nextCouplingGroupId; s.nextFaceEntryId = a.nextFaceEntryId; s.nextMatId = a.nextMatId; s.stepSurface = a.stepSurface; @@ -149,6 +160,7 @@ const createAnalysisActions: SliceCreator = (set) => ({ s.pickTargetGroupId = null; s.pickTieSide = "a"; s.tieDraft = { a: [], b: [] }; + s.couplingDraft = null; s.hasStarted = true; s.fitViewTrigger++; }), @@ -162,12 +174,14 @@ const createAnalysisActions: SliceCreator = (set) => ({ s.bcGroups = []; s.loadGroups = []; s.tieGroups = []; + s.couplingGroups = []; s.constraints = []; s.loads = []; s.surfaceLoads = []; s.nextBcGroupId = 1; s.nextLoadGroupId = 1; s.nextTieGroupId = 1; + s.nextCouplingGroupId = 1; s.nextFaceEntryId = 1; s.modelName = ""; s.result = null; @@ -191,6 +205,7 @@ const createAnalysisActions: SliceCreator = (set) => ({ s.pickTargetGroupId = null; s.pickTieSide = "a"; s.tieDraft = { a: [], b: [] }; + s.couplingDraft = null; s.mode = "geometry"; s.stepImportError = null; s.isRunning = false; diff --git a/web/src/wasm/pkg/kofem_wasm.d.ts b/web/src/wasm/pkg/kofem_wasm.d.ts index 8eb66fe..bdf6f7c 100644 --- a/web/src/wasm/pkg/kofem_wasm.d.ts +++ b/web/src/wasm/pkg/kofem_wasm.d.ts @@ -118,12 +118,16 @@ export interface KofemModule { * `attributes` is one 1-based solid-material index per TET, into * the `mat_json.solid` array; omitted means every tet uses the * first material. - * coupling: { ref, offsets, solid, mpc?, relaxation? } CSR-style: reference - * node ref[k] ties to solid[offsets[k]..offsets[k+1]). Optional - * mpc[k] = 1 selects the relaxed shell-to-solid MPC coupling (rigid - * translation tie + relaxation-scaled rotation) for that reference, + * coupling: { ref, offsets, solid, mpc?, dof_mask?, relaxation? } CSR-style: + * reference node ref[k] ties to solid[offsets[k]..offsets[k+1]). + * Optional mpc[k] = 1 selects the relaxed shell-to-solid MPC + * coupling (rigid translation tie + relaxation-scaled rotation), + * 2 the kinematic RBE2 coupling (the coupled nodes follow the + * reference point rigidly, leaving the point independent), and * 0/absent keeps the distributing RBE3 coupling; relaxation is the - * shared ψ ∈ [0.5,1] used by the MPC couplings. + * shared ψ ∈ [0.5,1] used by the MPC couplings. dof_mask[k] selects + * which of a kinematic coupling's six DOFs are tied (bits 0..5, + * all six when absent). * bcs: { fixed_dofs, load_dofs, load_vals } (DOF = 6·node+component) * mat_json: { solid, shell:{young_modulus,poisson_ratio} } — `solid` is * either one material object (every tet uses it) or an ARRAY of @@ -142,6 +146,7 @@ export interface KofemModule { offsets: Int32Array solid: Int32Array mpc?: Int32Array + dof_mask?: Int32Array relaxation?: number }, bcs: { fixed_dofs: Int32Array; load_dofs: Int32Array; load_vals: Float64Array }, diff --git a/web/src/workers/solver.worker.ts b/web/src/workers/solver.worker.ts index 03a0504..6202ca8 100644 --- a/web/src/workers/solver.worker.ts +++ b/web/src/workers/solver.worker.ts @@ -26,10 +26,18 @@ import { dropCouplingsOnFixedNodes, shellNodeLocator, isShellPoolIndex, + concatCouplings, + couplingMpcCodes, + couplingDofMasks, type CoupledModel, type ShellizeMesh, type ShellExtraction, } from "../lib/shellize.js"; +import { + buildReferenceCouplings, + referencePointIds, + type CouplingDefinition, +} from "../lib/coupling.js"; import { detectShellBodies } from "../lib/thinBodies.js"; let Module: KofemModule | null = null; @@ -192,6 +200,11 @@ interface SolvePayload { // parts that touch without a shared face are joined (#359). Absent/empty = // no ties, and the mesh reaches the solver untouched. tieGroups?: TieDefinition[]; + // Surface-to-point couplings (KOF-208): each idealises a picked surface to one + // reference point, either distributing (RBE3) or kinematic (RBE2). A model + // carrying any of these solves through the COUPLED assembler, whether or not + // it has shells. Absent/empty = no couplings, and the routing is unchanged. + couplings?: CouplingDefinition[]; // Surface mesh + per-triangle CAD face id (from meshing / the analysis file): // needed to detect thin-walled bodies and idealise them as shells (auto-shell). surfaceTriangles?: [number, number, number][] | null; @@ -1168,15 +1181,24 @@ function coupledMaterials( // clamped in rotation. `isShell` reports whether a pool node carries shell // (6-DOF) stiffness — the auto-shell and mixed paths supply it differently, but // the rule is the same. +// +// `isRefPoint` marks the pool nodes that are a coupling's REFERENCE POINT. Those +// carry six real DOFs (shell_core gives a coupling reference its rotations), so +// an Rx/Ry/Rz constraint on one is a genuine rotational restraint and is passed +// through instead of dropped — clamping a kinematic reference point is how a +// bolted connection is stated. Everywhere else a rotational constraint is still +// dropped: the shell nodes take their rotational clamp from the all-three- +// translations rule below, and a solid node has no rotational DOF to restrain. function coupledFixedDofs( constraints: Constraint[], poolOf: (nodeId: number) => number, isShell: (poolIndex: number) => boolean, + isRefPoint: (poolIndex: number) => boolean, ): number[] { const fixedByPool = new Map>(); for (const c of constraints) { - if (c.dof > 2) continue; const pi = poolOf(c.nodeId); + if (c.dof > 2 && !isRefPoint(pi)) continue; let dofs = fixedByPool.get(pi); if (!dofs) { dofs = new Set(); @@ -1187,23 +1209,43 @@ function coupledFixedDofs( const fixed_dofs: number[] = []; for (const [pi, dofs] of fixedByPool) { for (const d of dofs) fixed_dofs.push(6 * pi + d); - if (isShell(pi) && dofs.has(0) && dofs.has(1) && dofs.has(2)) + // A shell node clamped in all three translations is clamped, not hinged. + // A reference point is NOT given that treatment: its rotations are the DOFs + // the coupled surface's rigid-body motion rides on, so fixing them because + // the translations were fixed would silently turn a pinned point into a + // built-in one. Check Rx/Ry/Rz to clamp a reference point. + if ( + isShell(pi) && + !isRefPoint(pi) && + dofs.has(0) && + dofs.has(1) && + dofs.has(2) + ) for (const d of [3, 4, 5]) fixed_dofs.push(6 * pi + d); } return fixed_dofs; } // Point + surface loads → equivalent nodal forces on the pool. +// +// A moment (DOF 3..5) is applied directly where the node has a rotational DOF to +// receive it — a coupling REFERENCE POINT. That is the couple a surface-to-point +// coupling exists to carry: the point takes M, the coupling spreads it over the +// gripped surface as the statically equivalent traction, and no lever arm has to +// be invented. Elsewhere a moment still arrives pre-converted to a ring of nodal +// forces (momentToNodalForces), so it is skipped here as it always was. function coupledLoads( loads: Load[], surfaceLoads: SurfaceLoad[] | undefined, poolOf: (nodeId: number) => number, + isRefPoint: (poolIndex: number) => boolean, ): { load_dofs: number[]; load_vals: number[] } { const load_dofs: number[] = []; const load_vals: number[] = []; for (const l of loads) { - if (l.dof > 2) continue; - load_dofs.push(6 * poolOf(l.nodeId) + l.dof); + const pi = poolOf(l.nodeId); + if (l.dof > 2 && !isRefPoint(pi)) continue; + load_dofs.push(6 * pi + l.dof); load_vals.push(l.value); } for (const sl of surfaceLoads ?? []) { @@ -1234,7 +1276,8 @@ function mapCoupledDisplacements( ): Float64Array { const displacements = new Float64Array(3 * nodes.length); for (let i = 0; i < nodes.length; i++) { - const sp = model.solidPool.get(i); + const rp = model.refPool.get(i); + const sp = rp !== undefined ? rp : model.solidPool.get(i); const pi = sp !== undefined ? sp @@ -1356,14 +1399,21 @@ function tryCoupledSolve( // gapped pin/hole interface is a proper force-and-moment tie instead of a // sparse near-hinge. const wallTets = shellWallTets(mesh, shells); + const couplings = payload.couplings ?? []; + const refIds = referencePointIds(couplings); const model = buildCoupledModel(mesh, shells, wallTets, { ties: coupledTies(payload.tieGroups, elements, vid), + referencePoints: [...refIds].map((nodeId) => + vid(nodeId, "coupling reference point"), + ), }); - if (model.coupling.ref.length === 0) return null; // shell doesn't couple to the solid + if (model.coupling.ref.length === 0 && couplings.length === 0) return null; // shell doesn't couple to the solid const nearestShell = shellNodeLocator(model); const poolOf = (nodeId: number): number => { const vi = vid(nodeId, "coupled bc"); + const rp = model.refPool.get(vi); + if (rp !== undefined) return rp; const sp = model.solidPool.get(vi); if (sp !== undefined) return sp; return nearestShell([ @@ -1372,15 +1422,31 @@ function tryCoupledSolve( mesh.V[3 * vi + 2], ]); }; + const refPoolIndices = new Set(model.refPool.values()); + const isRefPoint = (pi: number) => refPoolIndices.has(pi); - const fixed_dofs = coupledFixedDofs(constraints, poolOf, (pi) => - isShellPoolIndex(model, pi), + const fixed_dofs = coupledFixedDofs( + constraints, + poolOf, + (pi) => isShellPoolIndex(model, pi), + isRefPoint, ); // A clamped shell rim can sit next to the retained base solid; the proximity // detector would otherwise couple the very nodes the user fixed (engine refuses // a fixed coupling-dependent node, #377). The BC wins. - const coupling = dropCouplingsOnFixedNodes(model.coupling, fixed_dofs); - const { load_dofs, load_vals } = coupledLoads(loads, surfaceLoads, poolOf); + // The declared surface-to-point couplings ride on the same pool mapping as the + // BCs: a coupled node whose thin wall was idealised away resolves to the + // mid-surface node that replaced it, which is where its stiffness now lives. + const coupling = concatCouplings( + dropCouplingsOnFixedNodes(model.coupling, fixed_dofs), + buildReferenceCouplings(couplings, (nodeId) => poolOf(nodeId)), + ); + const { load_dofs, load_vals } = coupledLoads( + loads, + surfaceLoads, + poolOf, + isRefPoint, + ); // Solid bodies that actually contribute solid tets to the pool: the other // bodies plus the shelled body when its thick base survived (only its thin walls @@ -1413,7 +1479,8 @@ function tryCoupledSolve( ref: Int32Array.from(coupling.ref), offsets: Int32Array.from(coupling.offsets), solid: Int32Array.from(coupling.solid), - mpc: Int32Array.from(coupling.mpc ?? coupling.ref.map(() => 0)), + mpc: Int32Array.from(couplingMpcCodes(coupling)), + dof_mask: Int32Array.from(couplingDofMasks(coupling)), relaxation: SHELL_SOLID_MPC_RELAXATION, }, { @@ -1589,7 +1656,16 @@ function tryPureShellSolve( // Reuse the coupled BC/load mapping (nearest-node lumping, 6·node+dof), then // fold it into the shell solver's fixed_vertices/fixed_dofs/point_loads form. - const flatFixed = coupledFixedDofs(constraints, shellOf, () => true); + // Every node here is a shell node and none is a coupling reference point: an + // all-shell model that declares a coupling is refused in handleSolve, because + // only the coupled assembler applies one. + const noReferencePoints: (poolIndex: number) => boolean = () => false; + const flatFixed = coupledFixedDofs( + constraints, + shellOf, + () => true, + noReferencePoints, + ); const dofsByVertex = new Map>(); for (const d of flatFixed) getOrInitDofs(dofsByVertex, Math.floor(d / 6)).add(d % 6); @@ -1600,7 +1676,12 @@ function tryPureShellSolve( else fixed_dofs.push({ vertex, dofs: [...dofSet].sort((a, b) => a - b) }); } - const { load_dofs, load_vals } = coupledLoads(loads, surfaceLoads, shellOf); + const { load_dofs, load_vals } = coupledLoads( + loads, + surfaceLoads, + shellOf, + noReferencePoints, + ); const forceByVertex = new Map(); for (let k = 0; k < load_dofs.length; k++) { const vertex = Math.floor(load_dofs[k] / 6); @@ -1747,6 +1828,27 @@ function mixedCoupledMaterials( const mat = materialOf(el, "shell"); shellUsed.set(mat.id, mat); } + // A model with NO shell element takes this path when it carries a + // surface-to-point coupling: the coupled assembler is what applies an RBE2/RBE3 + // constraint, and it accepts zero triangles. `mat.shell` is still read by + // solve_coupled, so it gets the first solid material — with no facet to + // assemble it can only be unused, and inventing an arbitrary modulus would + // print a stiffness that is not in the model anywhere. + if (shellUsed.size === 0 && shellElements.length === 0) + return { + mat: { + solid: solidOrder.map((mat) => ({ + young_modulus: mat.young, + poisson_ratio: mat.poisson, + })), + shell: { + young_modulus: solidOrder[0].young, + poisson_ratio: solidOrder[0].poisson, + }, + }, + solidAttributes, + solidMaterialNames: solidOrder.map((mat) => mat.name), + }; if (shellUsed.size === 0) throw new Error("mixed solve: the model has no shell element"); if (shellUsed.size > 1) { @@ -1802,8 +1904,12 @@ function resolveMixedThicknesses( // auto-shell path uses, only with the shells given explicitly instead of // idealised from thin solid walls. Constraints/loads map onto the 6-DOF pool the // same way the auto-shell path does (coupledFixedDofs/coupledLoads). elementOrder -// and tie connections do not apply (the coupled assembler is linear tets + DKT -// facets, whose interfaces are joined by RBE3 coupling) and are ignored. +// does not apply (the coupled assembler is linear tets + DKT facets). +// +// An ALL-SOLID model comes here too when it carries a surface-to-point coupling: +// an RBE2/RBE3 constraint only exists in this assembler, and it is happy with +// zero shell triangles. Nothing else about the path changes — the shell arrays +// are simply empty. function handleMixedSolve(id: number, payload: SolvePayload) { const { nodes, @@ -1853,12 +1959,20 @@ function handleMixedSolve(id: number, payload: SolvePayload) { } const thicknesses = resolveMixedThicknesses(shellElements, properties); + // Reference points belong to no element, so they enter the pool explicitly. + const couplings = payload.couplings ?? []; + const refIds = referencePointIds(couplings); const model = buildExplicitCoupledModel( verts, solidTets, shellTris, thicknesses, - { ties: coupledTies(payload.tieGroups, elements, vid) }, + { + ties: coupledTies(payload.tieGroups, elements, vid), + referencePoints: [...refIds].map((nodeId) => + vid(nodeId, "coupling reference point"), + ), + }, ); const poolOf = (nodeId: number): number => { @@ -1870,14 +1984,35 @@ function handleMixedSolve(id: number, payload: SolvePayload) { ); return pi; }; - const fixed_dofs = coupledFixedDofs(constraints, poolOf, (pi) => - model.shellPoolIndex.has(pi), + const refPoolIndices = new Set([...refIds].map((nodeId) => poolOf(nodeId))); + const isRefPoint = (pi: number) => refPoolIndices.has(pi); + const fixed_dofs = coupledFixedDofs( + constraints, + poolOf, + (pi) => model.shellPoolIndex.has(pi), + isRefPoint, ); // A clamped shell node that also sits within coupling range of the solid would // be both fixed and a distributing-coupling dependent — the engine refuses that // (#377). The BC wins; drop the coupling on those nodes. - const coupling = dropCouplingsOnFixedNodes(model.coupling, fixed_dofs); - const { load_dofs, load_vals } = coupledLoads(loads, surfaceLoads, poolOf); + const coupling = concatCouplings( + dropCouplingsOnFixedNodes(model.coupling, fixed_dofs), + buildReferenceCouplings(couplings, (nodeId, context) => { + const pi = model.poolOfVertex.get(vid(nodeId, context)); + if (pi === undefined) + throw new Error( + `${context}: node ${nodeId} is not part of the solved model — re-pick the ` + + "coupled surface, or delete the coupling.", + ); + return pi; + }), + ); + const { load_dofs, load_vals } = coupledLoads( + loads, + surfaceLoads, + poolOf, + isRefPoint, + ); const { mat, solidAttributes, solidMaterialNames } = mixedCoupledMaterials( materials, properties, @@ -1887,7 +2022,7 @@ function handleMixedSolve(id: number, payload: SolvePayload) { self.postMessage({ id, - log: `[mixed] ${solidElements.length} solid tets (${solidMaterialNames.join(", ")}), ${shellElements.length} shell facets → ${model.pool.length / 3} pool nodes, ${coupling.ref.length} couplings…`, + log: `[mixed] ${solidElements.length} solid tets (${solidMaterialNames.join(", ")}), ${shellElements.length} shell facets → ${model.pool.length / 3} pool nodes, ${coupling.ref.length} couplings (${couplings.length} declared)…`, }); const result = m().solve_coupled( @@ -1902,7 +2037,8 @@ function handleMixedSolve(id: number, payload: SolvePayload) { ref: Int32Array.from(coupling.ref), offsets: Int32Array.from(coupling.offsets), solid: Int32Array.from(coupling.solid), - mpc: Int32Array.from(coupling.mpc ?? coupling.ref.map(() => 0)), + mpc: Int32Array.from(couplingMpcCodes(coupling)), + dof_mask: Int32Array.from(couplingDofMasks(coupling)), relaxation: SHELL_SOLID_MPC_RELAXATION, }, { @@ -1955,11 +2091,32 @@ function handleSolve(id: number, payload: SolvePayload) { // explicit here, unlike the auto-shell path below, which idealises thin SOLID // bodies itself). const nShells = payload.elements.filter((e) => e.type === "CTRIA3").length; + // eslint-disable-next-line kofem/no-silent-fallback -- `couplings` is optional in the solve message; a model with none is the ordinary case + const nCouplings = payload.couplings?.length ?? 0; if (nShells > 0) { + // A coupling is an RBE2/RBE3 constraint, and only the coupled assembler + // applies one — the pure Kirchhoff shell solver has no couplings, and the + // coupled assembler needs a solid domain to assemble. An all-shell model + // therefore cannot carry a coupling today; say so, rather than solve it + // as if the coupling were not there. + if (nShells === payload.elements.length && nCouplings > 0) + throw new Error( + `This model is all shell elements and declares ${nCouplings} surface-to-point ` + + "coupling(s), which the shell solver cannot apply — a coupling needs the " + + "coupled solid-shell assembler, and that needs at least one solid body. " + + "Delete the coupling, or keep one body solid.", + ); if (nShells === payload.elements.length) handleShellSolve(id, payload); else handleMixedSolve(id, payload); return; } + // An all-solid model that declares a surface-to-point coupling still needs the + // coupled assembler — solve_linear_elastic has no notion of an RBE2/RBE3 + // constraint, and the reference point is not even a node it could carry. + if (nCouplings > 0) { + handleMixedSolve(id, payload); + return; + } // A body marked "Shell" is idealised as shells and solved coupled to the solid // bodies — this converges where the all-solid solve of the thin part stalls diff --git a/web/tests/couplings.spec.ts b/web/tests/couplings.spec.ts new file mode 100644 index 0000000..9314197 --- /dev/null +++ b/web/tests/couplings.spec.ts @@ -0,0 +1,391 @@ +// SPDX-FileCopyrightText: 2026 Michael Kofler +// SPDX-License-Identifier: AGPL-3.0-or-later + +// Surface-to-point couplings in the Constraints panel (KOF-208): a coupling is +// a named model object with a picked surface, a kind, a DOF mask and a +// reference point, created, edited and deleted like a BC or a load. +// +// The reference point is a real node of the model — created with the coupling, +// removed with it — which is what lets a BC or a load act on it through the +// ordinary face-entry machinery. These tests pin that lifecycle, because a +// leftover reference point is a node in no element, and the all-solid solve +// would assemble it into a singular system. +// +// Face picking itself needs 3D clicks, so the picked faces are written into the +// same pick-session state the viewport writes (selectedFace); everything from +// the kind select onwards is driven through the real UI. + +import { test, expect } from "./coverage"; +import type { Page } from "@playwright/test"; +import { boxHexMesh } from "../../examples/validation/lib/mesh.mjs"; + +interface CouplingGroupState { + id: number; + name: string; + kind: "distributing" | "kinematic"; + dofs: number[]; + refNodeId: number; + point: [number, number, number]; + faces: { id: number; nodeIds: number[] }[]; +} + +type Store = { + getState(): { + couplingGroups: CouplingGroupState[]; + bcGroups: { id: number; faces: { nodeIds: number[] }[] }[]; + loadGroups: { + id: number; + faces: { nodeIds: number[]; geometry?: string }[]; + }[]; + nodes: { id: number; x: number; y: number; z: number }[]; + couplingDraft: { + point: [number, number, number]; + nodeIds: number[]; + } | null; + selectedFace: { nodeIds: number[] } | null; + updateCouplingGroup( + id: number, + kind: string, + dofs: number[], + point: [number, number, number], + ): void; + setSelectedFace(face: { + nodeIds: number[]; + label: string; + axis: "X" | "Y" | "Z"; + isMax: boolean; + }): void; + }; + setState(s: object): void; +}; + +// One hex box, 10 × 4 × 4 mm. Its x = 10 face is the surface a coupling grips. +function box() { + const mesh = boxHexMesh(10, 4, 4, 5, 2, 2); + const nodes = mesh.vertices.map((v: number[], i: number) => ({ + id: i, + x: v[0], + y: v[1], + z: v[2], + })); + const elements = mesh.hexahedra.map((h: number[], i: number) => ({ + id: i, + type: "CHEXA", + nodeIds: [...h], + propertyId: 1, + })); + const endFace = nodes + .filter((node) => Math.abs(node.x - 10) < 1e-9) + .map((node) => node.id); + return { nodes, elements, endFace }; +} + +async function openConstraintsMode(page: Page) { + const model = box(); + await page.goto("/app/"); + await page.waitForFunction(() => + Boolean((window as unknown as { __kofemStore?: unknown }).__kofemStore), + ); + await page.evaluate((injected) => { + (window as unknown as { __kofemStore: Store }).__kofemStore.setState({ + nodes: injected.nodes, + elements: injected.elements, + properties: [{ id: 1, materialId: 1 }], + modelName: "Box", + hasStarted: true, + viewRepr: "surface", + mode: "constraints", + }); + }, model); + return model; +} + +// Write a picked face into the live pick session, the way a viewport click does. +async function pickFace(page: Page, nodeIds: number[]) { + await page.evaluate((ids) => { + (window as unknown as { __kofemStore: Store }).__kofemStore + .getState() + .setSelectedFace({ + nodeIds: ids, + label: `Face (${ids.length} nodes)`, + axis: "X", + isMax: true, + }); + }, nodeIds); +} + +const state = (page: Page) => + page.evaluate(() => + (window as unknown as { __kofemStore: Store }).__kofemStore.getState(), + ); + +test("a coupling idealises a picked surface to a reference point", async ({ + page, +}) => { + const model = await openConstraintsMode(page); + + await page.getByTestId("add-coupling").click(); + await pickFace(page, model.endFace); + + // The point defaults to the picked selection's centre — the x = 10 face of a + // 10 × 4 × 4 box, so (10, 2, 2). + await expect(page.getByTestId("coupling-point-x")).toHaveValue("10"); + await expect(page.getByTestId("coupling-point-y")).toHaveValue("2"); + await expect(page.getByTestId("coupling-point-z")).toHaveValue("2"); + + // Kinematic with the rotations dropped: a spherical joint at the point. + await page.getByTestId("coupling-kind").selectOption("kinematic"); + for (const dof of ["Rx", "Ry", "Rz"]) + await page.getByRole("checkbox", { name: dof, exact: true }).uncheck(); + await page.getByTestId("apply-coupling").click(); + + const { couplingGroups, nodes } = await state(page); + expect(couplingGroups).toHaveLength(1); + expect(couplingGroups[0].name).toBe("Coupling1"); + expect(couplingGroups[0].kind).toBe("kinematic"); + expect(couplingGroups[0].dofs).toEqual([0, 1, 2]); + expect(couplingGroups[0].point).toEqual([10, 2, 2]); + expect(couplingGroups[0].faces[0].nodeIds).toEqual(model.endFace); + + // The reference point is a real node, added to the model. + const refNode = nodes.find((n) => n.id === couplingGroups[0].refNodeId); + expect(refNode).toBeTruthy(); + expect([refNode?.x, refNode?.y, refNode?.z]).toEqual([10, 2, 2]); + expect(nodes).toHaveLength(model.nodes.length + 1); + + await expect(page.getByText("Coupling1")).toBeVisible(); + await expect( + page.getByText(/kinematic · Ux, Uy, Uz · 9 nodes/), + ).toBeVisible(); +}); + +test("a coupling's kind, DOFs and point are editable, and deleting it removes the reference point", async ({ + page, +}) => { + const model = await openConstraintsMode(page); + + await page.getByTestId("add-coupling").click(); + await pickFace(page, model.endFace); + await page.getByTestId("apply-coupling").click(); + + // ✎ opens the inline editor. Switching to distributing drops the DOF mask — + // a distributing coupling ties all six of its point's DOFs by construction, + // so the checkboxes are not offered and the stored mask is the full set. + await page.getByTitle("Edit coupling").click(); + const form = page.getByTestId("coupling-edit-form"); + await expect(form).toBeVisible(); + await form.getByTestId("coupling-kind").selectOption("distributing"); + await expect( + form.getByRole("checkbox", { name: "Rx", exact: true }), + ).toHaveCount(0); + await form.getByTestId("coupling-point-x").fill("14"); + await form.getByRole("button", { name: "Save" }).click(); + await expect(form).not.toBeVisible(); + + const edited = await state(page); + expect(edited.couplingGroups[0].kind).toBe("distributing"); + expect(edited.couplingGroups[0].dofs).toEqual([0, 1, 2, 3, 4, 5]); + expect(edited.couplingGroups[0].point).toEqual([14, 2, 2]); + // The reference point NODE moved with it — the point is the node. + const moved = edited.nodes.find( + (n) => n.id === edited.couplingGroups[0].refNodeId, + ); + expect(moved?.x).toBe(14); + + await page.getByTitle("Delete coupling").click(); + const after = await state(page); + expect(after.couplingGroups).toHaveLength(0); + expect(after.nodes).toHaveLength(model.nodes.length); +}); + +test("a BC can act on a coupling's reference point, and goes when the coupling does", async ({ + page, +}) => { + const model = await openConstraintsMode(page); + + await page.getByTestId("add-coupling").click(); + await pickFace(page, model.endFace); + await page.getByTestId("coupling-kind").selectOption("kinematic"); + await page.getByTestId("apply-coupling").click(); + const { couplingGroups } = await state(page); + const refNodeId = couplingGroups[0].refNodeId; + + // The reference point is picked in the viewport like any other selection — + // its marker carries the click (useReferencePointPick), which lands in the + // same session state a face pick writes. Clamping a kinematic point is how a + // bolted hole is restrained without fixing every node of its bore. + await page.getByRole("button", { name: "+ Add BC" }).click(); + await page.getByTestId("pick-geometry-point").click(); + await pickFace(page, [refNodeId]); + + // Picking a reference point exposes the rotational DOFs, because it has them + // — an ordinary solid node does not. + await expect( + page.getByRole("checkbox", { name: "Rx", exact: true }), + ).toBeVisible(); + await page.getByRole("button", { name: "Apply BC" }).click(); + + const withBc = await state(page); + expect(withBc.bcGroups).toHaveLength(1); + expect(withBc.bcGroups[0].faces[0].nodeIds).toEqual([refNodeId]); + + // Deleting the coupling takes the reference point with it, so the BC that + // named it must go too rather than constrain a node that no longer exists. + await page.getByTitle("Delete coupling").click(); + const after = await state(page); + expect(after.couplingGroups).toHaveLength(0); + expect(after.bcGroups).toHaveLength(0); + expect(after.nodes).toHaveLength(model.nodes.length); +}); + +test("the reference point being placed is previewed in the viewport", async ({ + page, +}) => { + const model = await openConstraintsMode(page); + + // Nothing is being placed yet, so there is nothing to preview. + await page.getByTestId("add-coupling").click(); + expect((await state(page)).couplingDraft).toBeNull(); + + // Picking the surface places the point at its centre, and the preview follows + // — the marker and the lines to the gripped nodes are what make it possible to + // judge the position before applying. + await pickFace(page, model.endFace); + const picked = await state(page); + expect(picked.couplingDraft?.point).toEqual([10, 2, 2]); + expect(picked.couplingDraft?.nodeIds).toEqual(model.endFace); + + // Typing coordinates moves the preview with them. + await page.getByTestId("coupling-point-x").fill("14"); + await expect + .poll(async () => (await state(page)).couplingDraft?.point[0]) + .toBe(14); + + // Half-typed input is not a position: the preview drops rather than park the + // marker at a number the user did not mean. + await page.getByTestId("coupling-point-x").fill(""); + await expect.poll(async () => (await state(page)).couplingDraft).toBeNull(); + + // Cancelling the pick session abandons the coupling, and the preview with it. + await page.getByTestId("coupling-point-x").fill("14"); + await expect + .poll(async () => (await state(page)).couplingDraft) + .not.toBeNull(); + await page.getByTitle("Cancel").click(); + expect((await state(page)).couplingDraft).toBeNull(); +}); + +test("a load can be applied at a single picked node", async ({ page }) => { + const model = await openConstraintsMode(page); + + await page.getByRole("button", { name: "+ Add Load" }).click(); + // A solid mesh offers Face and Point, but not Edge — a closed solid boundary + // has no boundary polyline to walk. + await expect(page.getByTestId("pick-geometry-face")).toBeVisible(); + await expect(page.getByTestId("pick-geometry-point")).toBeVisible(); + await expect(page.getByTestId("pick-geometry-edge")).toHaveCount(0); + + await page.getByTestId("pick-geometry-point").click(); + await pickFace(page, [model.endFace[0]]); + await expect( + page.getByText("applied at the picked node as a concentrated force"), + ).toBeVisible(); + await page.getByRole("button", { name: "Apply Load" }).click(); + + // The entry records that it was a POINT, which is what tells the solver to + // apply the force at the node instead of integrating it over a surface. + const { loadGroups } = await state(page); + expect(loadGroups).toHaveLength(1); + expect(loadGroups[0].faces[0].nodeIds).toEqual([model.endFace[0]]); + expect(loadGroups[0].faces[0].geometry).toBe("point"); +}); + +test("a pressure or a moment on a single node of a solid mesh is refused", async ({ + page, +}) => { + const model = await openConstraintsMode(page); + + await page.getByRole("button", { name: "+ Add Load" }).click(); + await page.getByTestId("pick-geometry-point").click(); + await pickFace(page, [model.endFace[0]]); + + // A pressure is force per unit area, and a node has none. + await page.getByRole("combobox").first().selectOption("pressure"); + await page.getByRole("button", { name: "Apply Load" }).click(); + await expect(page.getByTestId("constraints-error")).toContainText( + "cannot act on a single node", + ); + + // A moment needs a rotational DOF, which a solid mesh node does not have — + // the message points at the reference point, which does. + await page.getByRole("combobox").first().selectOption("moment"); + await page.getByRole("button", { name: "Apply Load" }).click(); + await expect(page.getByTestId("constraints-error")).toContainText( + "reference point", + ); + + expect((await state(page)).loadGroups).toHaveLength(0); +}); + +test("a reference point is picked by clicking its marker in the viewport", async ({ + page, +}) => { + const model = await openConstraintsMode(page); + + await page.getByTestId("add-coupling").click(); + await pickFace(page, model.endFace); + // Off the mesh grid on purpose: the face centre (10, 2, 2) is itself a node of + // this coarse box, and two exactly coincident candidates would not test which + // one the click resolves to. + await page.getByTestId("coupling-point-y").fill("1"); + await page.getByTestId("coupling-point-z").fill("1"); + await page.getByTestId("apply-coupling").click(); + const refNodeId = (await state(page)).couplingGroups[0].refNodeId; + + // Real viewport clicks, not injected state: the reference point belongs to no + // surface, so the ray always reports the mesh behind it — this is what proves + // the point still wins when it is what was clicked. + await page.getByRole("button", { name: "+ Add Load" }).click(); + await page.getByTestId("pick-geometry-point").click(); + + const canvas = page.locator("canvas").first(); + const box = await canvas.boundingBox(); + if (!box) throw new Error("viewport canvas has no bounding box"); + const hitsReferencePoint = () => + page.evaluate((id) => { + const selected = ( + window as unknown as { __kofemStore: Store } + ).__kofemStore.getState().selectedFace; + return selected?.nodeIds.length === 1 && selected.nodeIds[0] === id; + }, refNodeId); + + // Sweep for the marker — where it lands on screen depends on the camera. + let marker: [number, number] | null = null; + for (let fx = 0.55; fx <= 0.9 && !marker; fx += 0.04) { + for (let fy = 0.4; fy <= 0.85 && !marker; fy += 0.04) { + const spot: [number, number] = [ + box.x + box.width * fx, + box.y + box.height * fy, + ]; + await page.mouse.click(...spot); + if (await hitsReferencePoint()) marker = spot; + } + } + expect(marker).not.toBeNull(); + + // Start the session over and click the marker once, so the group carries the + // reference point and nothing the sweep collected on the way. + await page.getByTitle("Cancel").click(); + await page.getByRole("button", { name: "+ Add Load" }).click(); + await page.getByTestId("pick-geometry-point").click(); + if (!marker) throw new Error("no marker position"); + await page.mouse.click(...marker); + expect(await hitsReferencePoint()).toBe(true); + + await page.getByRole("button", { name: "Apply Load" }).click(); + const { loadGroups } = await state(page); + expect(loadGroups).toHaveLength(1); + expect(loadGroups[0].faces).toHaveLength(1); + expect(loadGroups[0].faces[0].nodeIds).toEqual([refNodeId]); + expect(loadGroups[0].faces[0].geometry).toBe("point"); +}); diff --git a/web/tests/edge-pick.spec.ts b/web/tests/edge-pick.spec.ts index 8783ad3..79ce762 100644 --- a/web/tests/edge-pick.spec.ts +++ b/web/tests/edge-pick.spec.ts @@ -13,9 +13,14 @@ import type { Page } from "@playwright/test"; type PickState = { selectedFace: { nodeIds: number[]; label: string } | null; - pickGeometry: "face" | "edge"; + pickGeometry: "face" | "edge" | "point"; bcGroups: { name: string; faces: { label: string; nodeIds: number[] }[] }[]; loadGroups: { name: string; faces: { label: string }[] }[]; + couplingGroups: { + name: string; + point: [number, number, number]; + faces: { label: string; nodeIds: number[] }[]; + }[]; nodes: unknown[]; }; @@ -43,6 +48,14 @@ function readPickState(page: Page) { name: g.name, faces: g.faces.map((f) => ({ label: f.label })), })), + couplingGroups: state.couplingGroups.map((g) => ({ + name: g.name, + point: g.point, + faces: g.faces.map((f) => ({ + label: f.label, + count: f.nodeIds.length, + })), + })), nodeCount: state.nodes.length, }; }); @@ -149,3 +162,44 @@ test("edge picking a shell rim creates an edge BC and an edge load", async ({ afterLoad.loadGroups[afterLoad.loadGroups.length - 1].faces[0]; expect(loadFace.label.startsWith("Edge")).toBe(true); }); + +test("a coupling can grip a picked shell edge", async ({ page }) => { + test.setTimeout(90_000); + + await page.goto("/app/?example=plate-with-hole-shell"); + await expect(page.locator("nav")).toBeVisible(); + await page.waitForFunction(() => + Boolean((window as unknown as { __kofem?: unknown }).__kofem), + ); + await expect + .poll(async () => (await readPickState(page)).nodeCount, { + timeout: 15_000, + }) + .toBeGreaterThan(0); + const totalNodes = (await readPickState(page)).nodeCount; + + await page.getByRole("button", { name: "Constraints" }).click(); + await page.getByTestId("add-coupling").click(); + + // A coupling grips a LINE as readily as a surface — the rim of a shell, a + // stiffener edge. On a flat sheet this is the only way to select one, because + // the whole sheet is a single face-pick region. + await page.getByTestId("pick-geometry-edge").click(); + await edgePickOnViewport(page); + + const picked = await readPickState(page); + expect(picked.pickGeometry).toBe("edge"); + expect(picked.selected?.label.startsWith("Edge")).toBe(true); + + // The reference point defaults to the centre of whatever was picked — for a + // rim, the centre of the ring. + await expect(page.getByTestId("coupling-point-x")).not.toHaveValue(""); + await page.getByTestId("apply-coupling").click(); + + const after = await readPickState(page); + expect(after.couplingGroups).toHaveLength(1); + const gripped = after.couplingGroups[0].faces[0]; + expect(gripped.label.startsWith("Edge")).toBe(true); + expect(gripped.count).toBeGreaterThan(2); + expect(gripped.count).toBeLessThan(totalNodes * 0.5); +}); diff --git a/web/tests/test_face_pick.mjs b/web/tests/test_face_pick.mjs index 68a3727..afe4540 100644 --- a/web/tests/test_face_pick.mjs +++ b/web/tests/test_face_pick.mjs @@ -16,6 +16,8 @@ import { buildBoundaryMeshTopo, mapTrianglesToCadFaces, pickFaceNodeIds, + pickPointNodeId, + nearestReferencePoint, toggleFaceSelection, } from "../src/lib/facePick.ts"; @@ -354,6 +356,80 @@ assert( triangles[holeOuter].every((v) => !holePick.has(v)), ); +// ── Point picking ───────────────────────────────────────────────────────────── +// +// A point pick selects ONE node: the corner of the clicked facet nearest the +// click. The candidates are restricted to that facet on purpose — a global +// nearest-node search would happily return a node on the far side of the model +// that projects closer in space, and the user would have selected something they +// cannot see. + +const pointTri = innerTriIndices[0]; +const [pa, pb, pc] = triangles[pointTri]; + +// Click exactly on one corner: that corner is the pick. +for (const corner of [pa, pb, pc]) { + const picked = pickPointNodeId(getPos(corner), pointTri, topoWithIds, getPos); + assert( + `a click on node ${corner} picks node ${corner}`, + setsEqual(picked, new Set([corner])), + ); +} + +// Click just off a corner, still inside the facet: the same corner wins. +const posA = getPos(pa); +const posB = getPos(pb); +const nearA = [ + posA[0] + 0.2 * (posB[0] - posA[0]), + posA[1] + 0.2 * (posB[1] - posA[1]), + posA[2] + 0.2 * (posB[2] - posA[2]), +]; +assert( + "a click 20 % along an edge picks the nearer end", + setsEqual( + pickPointNodeId(nearA, pointTri, topoWithIds, getPos), + new Set([pa]), + ), +); + +// Every pick is a single node, whatever the click — the invariant the load path +// depends on (a one-node selection spans no element face, so its load is +// applied at the node rather than integrated). +const centroid = [0, 1, 2].map( + (axis) => (getPos(pa)[axis] + getPos(pb)[axis] + getPos(pc)[axis]) / 3, +); +assert( + "a click at the facet centroid still picks exactly one node", + pickPointNodeId(centroid, pointTri, topoWithIds, getPos).size === 1, +); + +// ── Reference points compete with mesh nodes by distance ───────────────────── +// +// A reference point belongs to no surface, so the ray can only ever report the +// mesh behind it — a marker drawn on the face it couples is ALWAYS behind that +// face, whatever its render order. Which of the two was meant is therefore a +// distance question, and this is the comparison that answers it. + +const refA = { nodeId: 900, point: [0, 0, 0] }; +const refB = { nodeId: 901, point: [10, 0, 0] }; + +assert( + "the nearer of two reference points wins", + nearestReferencePoint([9, 0, 0], [refA, refB])?.nodeId === 901, +); +assert( + "…and the other way round", + nearestReferencePoint([1, 0, 0], [refA, refB])?.nodeId === 900, +); +assert( + "the squared distance is reported so the caller can compare against a node", + nearestReferencePoint([3, 4, 0], [refA])?.distSq === 25, +); +assert( + "a model with no couplings has no reference point to pick", + nearestReferencePoint([0, 0, 0], []) === null, +); + // ── Summary ─────────────────────────────────────────────────────────────────── console.log(`\n${passed + failed} tests: ${passed} passed, ${failed} failed`); diff --git a/web/tests/test_reference_point.mjs b/web/tests/test_reference_point.mjs new file mode 100644 index 0000000..6ab42a8 --- /dev/null +++ b/web/tests/test_reference_point.mjs @@ -0,0 +1,695 @@ +#!/usr/bin/env bun +// SPDX-FileCopyrightText: 2026 Michael Kofler +// SPDX-License-Identifier: AGPL-3.0-or-later + +// Surface-to-point coupling, from the model down to the engine (KOF-208, part 2). +// +// tests/test_coupling.mjs already pins the ENGINE's kinematic coupling against +// hand-built pool arrays. This one drives the same engine through the model-side +// pipeline the app actually uses — store-shaped nodes/elements, a coupling +// declared over a picked surface, buildExplicitCoupledModel + +// buildReferenceCouplings, then solve_coupled — so a break anywhere in that +// chain fails here rather than only in a browser. +// +// Checks: +// 1. Reference point placement. On a cylindrical bore, referencePointOptions +// offers the surface centre AND the two ends of the fitted axis (the "up +// and down centre" of KOF-208); on a flat face it offers the centre alone, +// because axis-end positions on a flat patch would be meaningless. +// 2. An all-solid model with a coupling reaches the engine at all: the model +// has no shell element, so this is the routing PR #417 left for part 2. +// 3. A force at the reference point of a KINEMATIC coupling produces the same +// tip deflection as the same force applied straight to the coupled face. +// 4. A MOMENT at the reference point bends the beam — the whole point of the +// feature, and impossible on a solid mesh without it (a solid node has no +// rotational DOF, and momentToNodalForces cannot build a couple from one +// node). rebuildLoads must turn it into a rotational DOF load. +// Both kinds of load on a point are pinned here, including that exactly one +// builder claims the face: rebuildLoads applies it as a nodal DOF load and +// rebuildSurfaceLoads leaves it alone, so it is neither dropped nor applied +// twice. +// 5. A kinematic reference point can be FIXED, which a distributing one +// cannot: clamping the point clamps the surface it grips. +// 6. Modelling errors are named, not left to the engine: a distributing +// coupling that grips fewer than 3 nodes, and two kinematic couplings +// fighting over the same node. +// +// Usage: bun tests/test_reference_point.mjs (from the web/ directory) + +import { readFileSync } from "fs"; +import { fileURLToPath } from "url"; +import { dirname, join } from "path"; +import { buildExplicitCoupledModel } from "../src/lib/shellize.ts"; +import { + buildReferenceCouplings, + referencePointOptions, +} from "../src/lib/coupling.ts"; +import { + rebuildLoads, + rebuildSurfaceLoads, +} from "../src/store/boundarySlice.ts"; + +const __dirname = dirname(fileURLToPath(import.meta.url)); +const wasmPkg = join(__dirname, "../src/wasm/pkg"); +const wasmBinary = readFileSync(join(wasmPkg, "kofem_wasm_emcc.wasm")).buffer; +const { default: createModule } = await import( + join(wasmPkg, "kofem_wasm_emcc.js") +); + +let failures = 0; +function check(name, cond, detail = "") { + if (cond) { + console.log(` [PASS] ${name}`); + } else { + failures++; + console.log(` [FAIL] ${name}${detail ? ` — ${detail}` : ""}`); + } +} + +const YOUNG = 210000; // MPa +const POISSON = 0.3; + +// ── Fixtures ────────────────────────────────────────────────────────────────── + +const KUHN = [ + [0, 1, 3, 7], + [0, 3, 2, 7], + [0, 2, 6, 7], + [0, 6, 4, 7], + [0, 4, 5, 7], + [0, 5, 1, 7], +]; + +// Store-shaped cantilever: a box [0,L] × [0,h] × [0,h] of Kuhn tets. +function beamModel(length, height, nx, nt = 2) { + const nodes = []; + const index = new Map(); + const dx = length / nx; + const dy = height / nt; + const at = (i, j, k) => { + const [x, y, z] = [i * dx, j * dy, k * dy]; + const key = `${x.toFixed(6)},${y.toFixed(6)},${z.toFixed(6)}`; + let id = index.get(key); + if (id === undefined) { + id = nodes.length; + index.set(key, id); + nodes.push({ id, x, y, z }); + } + return id; + }; + const elements = []; + for (let i = 0; i < nx; i++) + for (let j = 0; j < nt; j++) + for (let k = 0; k < nt; k++) { + const corner = [ + at(i, j, k), + at(i + 1, j, k), + at(i, j + 1, k), + at(i + 1, j + 1, k), + at(i, j, k + 1), + at(i + 1, j, k + 1), + at(i, j + 1, k + 1), + at(i + 1, j + 1, k + 1), + ]; + for (const tet of KUHN) + elements.push({ + id: elements.length, + type: "CTETRA", + nodeIds: [ + corner[tet[0]], + corner[tet[1]], + corner[tet[2]], + corner[tet[3]], + ], + propertyId: 1, + }); + } + return { nodes, elements }; +} + +// A hollow cylinder (annulus extruded along z) as store-shaped tets, so the +// bore is a real cylindrical surface the axis fit has to recognise. +function boreModel({ + radius = 10, + wall = 4, + height = 24, + nTheta = 16, + nz = 4, +}) { + const nodes = []; + const elements = []; + const id = (ring, theta, layer) => + layer * (2 * nTheta) + ring * nTheta + (theta % nTheta); + for (let layer = 0; layer <= nz; layer++) + for (let ring = 0; ring < 2; ring++) + for (let theta = 0; theta < nTheta; theta++) { + const angle = (2 * Math.PI * theta) / nTheta; + const ringRadius = ring === 0 ? radius : radius + wall; + nodes.push({ + id: nodes.length, + x: ringRadius * Math.cos(angle), + y: ringRadius * Math.sin(angle), + z: (height * layer) / nz, + }); + } + for (let layer = 0; layer < nz; layer++) + for (let theta = 0; theta < nTheta; theta++) { + const corner = [ + id(0, theta, layer), + id(0, theta + 1, layer), + id(1, theta, layer), + id(1, theta + 1, layer), + id(0, theta, layer + 1), + id(0, theta + 1, layer + 1), + id(1, theta, layer + 1), + id(1, theta + 1, layer + 1), + ]; + for (const tet of KUHN) + elements.push({ + id: elements.length, + type: "CTETRA", + nodeIds: [ + corner[tet[0]], + corner[tet[1]], + corner[tet[2]], + corner[tet[3]], + ], + propertyId: 1, + }); + } + return { nodes, elements }; +} + +// ── The solve, exactly as handleMixedSolve assembles it ────────────────────── + +const Module = await createModule({ wasmBinary }); + +function solve({ nodes, elements, couplings, constraints, loads }) { + const vertexIndex = new Map(nodes.map((node, i) => [node.id, i])); + const vid = (nodeId) => { + const index = vertexIndex.get(nodeId); + if (index === undefined) throw new Error(`unknown node ${nodeId}`); + return index; + }; + const verts = []; + for (const node of nodes) verts.push(node.x, node.y, node.z); + const solidTets = []; + for (const el of elements) + for (const nid of el.nodeIds) solidTets.push(vid(nid)); + + const refIds = [...new Set(couplings.map((c) => c.refNodeId))]; + const model = buildExplicitCoupledModel(verts, solidTets, [], [], { + referencePoints: refIds.map(vid), + }); + const poolOf = (nodeId) => { + const pi = model.poolOfVertex.get(vid(nodeId)); + if (pi === undefined) throw new Error(`node ${nodeId} is not in the pool`); + return pi; + }; + const refPool = new Set(refIds.map(poolOf)); + const coupling = buildReferenceCouplings(couplings, (nodeId) => + poolOf(nodeId), + ); + + const fixed = []; + for (const c of constraints) { + const pi = poolOf(c.nodeId); + if (c.dof > 2 && !refPool.has(pi)) continue; + fixed.push(6 * pi + c.dof); + } + const loadDofs = []; + const loadVals = []; + for (const l of loads) { + const pi = poolOf(l.nodeId); + if (l.dof > 2 && !refPool.has(pi)) continue; + loadDofs.push(6 * pi + l.dof); + loadVals.push(l.value); + } + + const result = Module.solve_coupled( + { + vertices: Float64Array.from(model.pool), + tets: Int32Array.from(model.tets), + triangles: Int32Array.from([]), + thicknesses: Float64Array.from([]), + }, + { + ref: Int32Array.from(coupling.ref), + offsets: Int32Array.from(coupling.offsets), + solid: Int32Array.from(coupling.solid), + mpc: Int32Array.from(coupling.mpc), + dof_mask: Int32Array.from(coupling.dofMask), + relaxation: 1.0, + }, + { + fixed_dofs: Int32Array.from(fixed), + load_dofs: Int32Array.from(loadDofs), + load_vals: Float64Array.from(loadVals), + }, + JSON.stringify({ + solid: { young_modulus: YOUNG, poisson_ratio: POISSON }, + shell: { young_modulus: YOUNG, poisson_ratio: POISSON }, + }), + ); + if ("error" in result) throw new Error(result.error); + return { result, poolOf }; +} + +// Largest |u| over a set of store node ids. +function maxDisplacement(result, poolOf, nodeIds) { + let best = 0; + for (const nodeId of nodeIds) { + const pi = poolOf(nodeId); + const mag = Math.hypot( + result.displacements[3 * pi], + result.displacements[3 * pi + 1], + result.displacements[3 * pi + 2], + ); + if (mag > best) best = mag; + } + return best; +} + +// One displacement component of one store node. +function displacement(result, poolOf, nodeId, dof) { + return result.displacements[3 * poolOf(nodeId) + dof]; +} + +// ── 1. Reference point placement ───────────────────────────────────────────── + +console.log("Reference point placement"); +{ + const bore = boreModel({}); + // The bore surface: the inner ring of nodes, which is exactly what a face + // pick on the hole would return. + const boreFace = { + nodeIds: bore.nodes + .filter((n) => Math.hypot(n.x, n.y) < 10 + 1e-6) + .map((n) => n.id), + }; + const options = referencePointOptions([boreFace], bore.nodes, bore.elements); + check( + "a cylindrical bore offers the centre plus both axis ends", + options.length === 3, + `got ${options.map((o) => o.label).join(", ")}`, + ); + const centre = options[0].point; + check( + "the selection centre is on the bore axis at mid-height", + Math.hypot(centre[0], centre[1]) < 1e-6 && Math.abs(centre[2] - 12) < 1e-6, + `got (${centre.join(", ")})`, + ); + const ends = options.slice(1).map((o) => o.point); + const zs = ends.map((p) => p[2]).sort((a, b) => a - b); + check( + "the axis ends sit at the two ends of the bore", + ends.every((p) => Math.hypot(p[0], p[1]) < 1e-6) && + Math.abs(zs[0]) < 1e-6 && + Math.abs(zs[1] - 24) < 1e-6, + `got z = ${zs.join(", ")}`, + ); + + // A flat face is not a cylinder — offering it "axis ends" would place the + // reference point somewhere with no meaning. + const beam = beamModel(100, 20, 5); + const endFace = { + nodeIds: beam.nodes.filter((n) => n.x > 100 - 1e-6).map((n) => n.id), + }; + const flatOptions = referencePointOptions( + [endFace], + beam.nodes, + beam.elements, + ); + check( + "a flat face offers only the selection centre", + flatOptions.length === 1 && flatOptions[0].label === "Selection centre", + `got ${flatOptions.map((o) => o.label).join(", ")}`, + ); +} + +// ── 2-4. An all-solid model with a coupling, loaded at its point ───────────── + +console.log("\nKinematic coupling on an all-solid cantilever"); +const LENGTH = 100; +const HEIGHT = 20; +const FORCE = -5000; // N, in −y + +const beam = beamModel(LENGTH, HEIGHT, 8); +const tipNodeIds = beam.nodes + .filter((node) => node.x > LENGTH - 1e-6) + .map((node) => node.id); +const rootNodeIds = beam.nodes + .filter((node) => node.x < 1e-6) + .map((node) => node.id); +// The reference point is a node of the model, as the store creates it. +const REF_ID = beam.nodes.length; +const refNode = { id: REF_ID, x: LENGTH, y: HEIGHT / 2, z: HEIGHT / 2 }; +const withRef = { ...beam, nodes: [...beam.nodes, refNode] }; + +const clamp = rootNodeIds.flatMap((nodeId) => + [0, 1, 2].map((dof) => ({ nodeId, dof })), +); +const coupling = { + name: "Coupling1", + kind: "kinematic", + dofs: [0, 1, 2, 3, 4, 5], + refNodeId: REF_ID, + faces: [{ nodeIds: tipNodeIds }], +}; + +// Centre of the tip face, on the beam's neutral axis. Deflection is read there +// rather than as a max over the face: under an end moment the face ROTATES, so +// a magnitude over its corners mixes bending deflection with that rotation. +const centreTipId = beam.nodes.find( + (node) => + node.x > LENGTH - 1e-6 && + Math.abs(node.y - HEIGHT / 2) < 1e-6 && + Math.abs(node.z - HEIGHT / 2) < 1e-6, +).id; + +// Tip deflection under the end force, on this mesh. Linear tets lock in +// bending, so the absolute value is well short of FL³/(3EI) at this element +// size — which is why the moment check below compares against THIS rather than +// against the analytic beam: the ratio of the two load cases isolates the +// coupling from the element's own stiffness error. +let forceDeflection = 0; + +{ + // The same total force, once through the reference point and once spread + // straight over the coupled face. + const throughPoint = solve({ + ...withRef, + couplings: [coupling], + constraints: clamp, + loads: [{ nodeId: REF_ID, dof: 1, value: FORCE }], + }); + const direct = solve({ + ...beam, + couplings: [], + constraints: clamp, + loads: tipNodeIds.map((nodeId) => ({ + nodeId, + dof: 1, + value: FORCE / tipNodeIds.length, + })), + }); + forceDeflection = Math.abs( + displacement(throughPoint.result, throughPoint.poolOf, centreTipId, 1), + ); + const uDirect = Math.abs( + displacement(direct.result, direct.poolOf, centreTipId, 1), + ); + check( + "an all-solid model with a coupling solves at all", + forceDeflection > 0, + `tip uy = ${forceDeflection}`, + ); + const relative = Math.abs(forceDeflection - uDirect) / uDirect; + check( + "a force at the reference point matches the force on the face", + relative < 0.05, + `${forceDeflection.toFixed(5)} vs ${uDirect.toFixed(5)} mm (${(100 * relative).toFixed(2)} %)`, + ); +} + +{ + // A MOMENT at the reference point. Only a rigid spider can produce one on a + // solid mesh: it is the θ_R × r term of u_i = u_R + θ_R × r_i that turns the + // couple into a self-equilibrated force pattern over the face. + // + // Compared against the force case on the same mesh, where the ratio is pure + // beam theory and independent of how stiff the discretisation is: + // δ_M/δ_F = [M·L²/(2EI)] / [F·L³/(3EI)] = 3M/(2FL) + const MOMENT = 2e5; // N·mm about z + const expectedRatio = (3 * MOMENT) / (2 * Math.abs(FORCE) * LENGTH); + const { result, poolOf } = solve({ + ...withRef, + couplings: [coupling], + constraints: clamp, + loads: [{ nodeId: REF_ID, dof: 5, value: MOMENT }], + }); + const tip = Math.abs(displacement(result, poolOf, centreTipId, 1)); + const ratio = tip / forceDeflection; + check( + "a moment at the reference point bends the beam", + tip > 0, + `tip uy = ${tip}`, + ); + const relative = Math.abs(ratio - expectedRatio) / expectedRatio; + check( + "the moment's tip deflection is 3M/(2FL) of the force case, as beam theory says", + relative < 0.05, + `ratio ${ratio.toFixed(4)} vs ${expectedRatio.toFixed(4)} (${(100 * relative).toFixed(2)} %)`, + ); + + // …and the model layer is what turns that moment group into a rotational DOF + // load. Without this, momentToNodalForces would see a one-node face, find + // every lever arm zero, and drop the load with a console warning. + const loads = rebuildLoads( + [ + { + id: 1, + name: "Load1", + dof: 5, + totalForce: MOMENT, + components: [0, 0, MOMENT], + kind: "moment", + faces: [ + { id: 1, label: "Coupling1 reference point", nodeIds: [REF_ID] }, + ], + }, + ], + withRef.nodes, + [ + { + ...coupling, + id: 1, + point: [refNode.x, refNode.y, refNode.z], + faces: [], + }, + ], + ); + check( + "a moment on a reference point becomes a rotational DOF load", + loads.length === 1 && loads[0].dof === 5 && loads[0].value === MOMENT, + JSON.stringify(loads), + ); + + // A FORCE on a reference point has the same problem from the other side: it + // would normally be integrated as a traction over the boundary faces of the + // selection, and a lone point spans none, so it too has to become a nodal DOF + // load rather than silently integrate to nothing. + const forceLoads = rebuildLoads( + [ + { + id: 2, + name: "Load2", + dof: 1, + totalForce: FORCE, + components: [0, FORCE, 0], + kind: "force", + faces: [ + { id: 2, label: "Coupling1 reference point", nodeIds: [REF_ID] }, + ], + }, + ], + withRef.nodes, + [ + { + ...coupling, + id: 1, + point: [refNode.x, refNode.y, refNode.z], + faces: [], + }, + ], + ); + check( + "a force on a reference point becomes a translational DOF load", + forceLoads.length === 1 && + forceLoads[0].dof === 1 && + forceLoads[0].value === FORCE, + JSON.stringify(forceLoads), + ); + + // …and the surface-load builder must NOT also claim that face, or the force + // would be applied twice. The two builders partition a group's faces by the + // same test, so this pins the split rather than trusting that the boundary + // matcher happens to find nothing on a one-node selection. + const pointForceGroup = { + id: 2, + name: "Load2", + dof: 1, + totalForce: FORCE, + components: [0, FORCE, 0], + kind: "force", + faces: [{ id: 2, label: "Coupling1 reference point", nodeIds: [REF_ID] }], + }; + const couplingGroups = [ + { ...coupling, id: 1, point: [refNode.x, refNode.y, refNode.z], faces: [] }, + ]; + check( + "a force on a reference point produces no surface load to double-apply it", + rebuildSurfaceLoads([pointForceGroup], withRef.elements, couplingGroups) + .length === 0, + ); + // A force on a real surface still takes the traction route, untouched. + check( + "a force on a picked surface still becomes a surface load", + rebuildSurfaceLoads( + [ + { + ...pointForceGroup, + faces: [{ id: 3, label: "Face 1", nodeIds: tipNodeIds }], + }, + ], + withRef.elements, + couplingGroups, + ).length === 1, + ); + + // A POINT pick — a single mesh node, marked `geometry: "point"` — goes the + // same way as a reference point, and for the same reason: one node spans no + // element face, so an integrated traction over it is a traction over nothing. + // The marker is explicit rather than inferred from the node count, so the two + // builders agree by construction instead of by coincidence. + const pointPickGroup = { + id: 3, + name: "Load3", + dof: 1, + totalForce: FORCE, + components: [0, FORCE, 0], + kind: "force", + faces: [ + { + id: 4, + label: "Node 1", + nodeIds: [centreTipId], + geometry: "point", + }, + ], + }; + const nodeLoads = rebuildLoads([pointPickGroup], withRef.nodes, []); + check( + "a force on a picked node becomes a nodal DOF load", + nodeLoads.length === 1 && + nodeLoads[0].nodeId === centreTipId && + nodeLoads[0].dof === 1 && + nodeLoads[0].value === FORCE, + JSON.stringify(nodeLoads), + ); + check( + "a force on a picked node produces no surface load to double-apply it", + rebuildSurfaceLoads([pointPickGroup], withRef.elements, []).length === 0, + ); +} + +{ + // Clamping the reference point clamps the surface it grips — the bolted-hole + // idealisation, and the thing a distributing coupling cannot do. + const { result, poolOf } = solve({ + ...withRef, + couplings: [coupling], + constraints: [ + ...clamp, + ...[0, 1, 2, 3, 4, 5].map((dof) => ({ nodeId: REF_ID, dof })), + ], + loads: rootNodeIds + .slice(0, 1) + .map((nodeId) => ({ nodeId, dof: 1, value: 0 })), + }); + const tip = maxDisplacement(result, poolOf, tipNodeIds); + check( + "a clamped kinematic reference point holds its coupled surface still", + tip < 1e-9, + `tip |u| = ${tip}`, + ); +} + +// ── 6. Modelling errors are named ──────────────────────────────────────────── + +console.log("\nModelling errors"); +{ + const poolOf = (nodeId) => nodeId; + let message = ""; + try { + buildReferenceCouplings( + [ + { + name: "Coupling1", + kind: "distributing", + dofs: [0, 1, 2, 3, 4, 5], + refNodeId: 99, + faces: [{ nodeIds: [1, 2] }], + }, + ], + poolOf, + ); + } catch (err) { + message = err.message; + } + check( + "a distributing coupling with fewer than 3 nodes is refused by name", + message.includes("Coupling1") && message.includes("at least 3"), + message, + ); + + message = ""; + try { + buildReferenceCouplings( + [ + { + name: "Bolt", + kind: "kinematic", + dofs: [0, 1, 2], + refNodeId: 90, + faces: [{ nodeIds: [1, 2, 3] }], + }, + { + name: "Bearing", + kind: "kinematic", + dofs: [0, 1, 2], + refNodeId: 91, + faces: [{ nodeIds: [3, 4, 5] }], + }, + ], + poolOf, + ); + } catch (err) { + message = err.message; + } + check( + "two kinematic couplings gripping the same node name both of them", + message.includes("Bolt") && message.includes("Bearing"), + message, + ); + + message = ""; + try { + buildReferenceCouplings( + [ + { + name: "Coupling1", + kind: "kinematic", + dofs: [], + refNodeId: 90, + faces: [{ nodeIds: [1, 2, 3] }], + }, + ], + poolOf, + ); + } catch (err) { + message = err.message; + } + check( + "a kinematic coupling that ties no DOF is refused", + message.includes("ties no DOF"), + message, + ); +} + +console.log( + failures === 0 + ? "\nAll reference-point checks passed" + : `\n${failures} check(s) FAILED`, +); +process.exit(failures === 0 ? 0 : 1);