From af7e282995b7774d0a9427fbf6d82a2ec3f5bdc8 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 19 Aug 2026 05:46:04 +0000 Subject: [PATCH] Support prescribed (non-zero) displacements on shell models The solid path gets inhomogeneous essential BCs from MFEM, which eliminates a seeded GridFunction in FormLinearSystem. The shell formulation had no counterpart: ShellInput.fixed_dofs are "constrained to zero" and neither solve_shell's nor solve_coupled's bcs object carried prescribed values at all. A displacement-driven shell model therefore refused to run (explicit CTRIA3) or, before it was made loud, silently pinned the value to zero and returned an all-zero field that looks like a converged answer. Thin bodies are marked as shells automatically at import, so a user meshing a thin-walled part lands on this path without ever choosing it. Engine: - ShellInput and CoupledInput gain prescribed_dofs (global DOF -> value), alongside the homogeneous fixed_dofs rather than replacing it. - apply_homogeneous_bc becomes apply_essential_bc: the known column contribution K[:,p]*g moves onto the right-hand side in a column pass, while the free rows still hold their original entries, before the row/column is cleared and F[p] = g. An empty value vector reproduces the old homogeneous behaviour exactly, so undriven models are untouched. - The coupling reduction carries prescribed values into the reduced system; an independent DOF maps one-to-one onto its reduced column, so the value needs no transform. - Plumbed through both Embind entry points: solve_shell takes bcs_json.prescribed_dofs [{vertex, dof, value}] (components 0..5, so a shell rotation is drivable too), solve_coupled takes bcs.prescribed_dofs and bcs.prescribed_vals. A prescribed DOF on a distributing-coupling reference node is refused the same way a fixed one already is (#377): the coupling, not the boundary condition, governs that node's motion either way. The worker drops an auto-detected coupling on a driven node exactly as it does on a clamped one, so the BC wins before the engine sees the conflict. Worker: both refusals replaced with the real mapping, splitting constraints into homogeneous and prescribed the way groupDirichlet already does for solids, across all three shell paths (explicit CTRIA3, pure auto-shell, coupled/mixed). A driven node is not given the all-translations-clamped -> rotations-clamped treatment, since a face pulled to a displacement is not thereby built in. A driven rotation on a node with no drivable rotational DOF now throws instead of being dropped by the rotational rule -- that silent drop is the same failure mode this change exists to remove. Verification: six native checks in engine/tests/shell_validation.cpp -- prescribed extension reaches delta, the driven edge still shows the correct Poisson contraction (driven, not clamped), the displacement drive agrees with the equivalent force drive (which is what proves the K[:,p]*g term), a coupled model driven through its RBE3 anchors translates rigidly, and both rejection paths throw. examples/validation/shell-prescribed- displacement.test.mjs mirrors the solid regression test through the WASM engine and is wired into bun run test. Fixes KOF-210 Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_012eF9TGvCgM9jiPfBEuqghx --- engine/cpp/shell_core.cpp | 99 +++++++--- engine/cpp/shell_core.h | 8 + engine/cpp/solve_coupled.cpp | 17 +- engine/cpp/solve_coupled.h | 5 +- engine/cpp/solve_shell.cpp | 25 ++- engine/cpp/solve_shell.h | 5 +- engine/tests/shell_validation.cpp | 176 ++++++++++++++++++ examples/validation/README.md | 23 +++ examples/validation/lib/solver.mjs | 59 +++++- .../shell-prescribed-displacement.test.mjs | 168 +++++++++++++++++ web/package.json | 2 +- web/src/wasm/pkg/kofem_wasm.d.ts | 21 ++- web/src/workers/solver.worker.ts | 157 ++++++++++++---- 13 files changed, 682 insertions(+), 83 deletions(-) create mode 100644 examples/validation/shell-prescribed-displacement.test.mjs diff --git a/engine/cpp/shell_core.cpp b/engine/cpp/shell_core.cpp index 56682611..a0da1ee7 100644 --- a/engine/cpp/shell_core.cpp +++ b/engine/cpp/shell_core.cpp @@ -555,20 +555,61 @@ std::array, 3> skew(const Vec3& a) { return {{{0.0, -a[2], a[1]}, {a[2], 0.0, -a[0]}, {-a[1], a[0], 0.0}}}; } -void apply_homogeneous_bc(Sparse& K, std::vector& F, const std::vector& fixed) { +// Eliminate the essential (Dirichlet) DOFs. `fixed` marks every constrained DOF +// and `values` carries what each is constrained TO — zero for a homogeneous +// clamp, g for a prescribed displacement. Standard inhomogeneous elimination: +// the known column contribution K[:,p]·g moves to the right-hand side before the +// row/column is cleared, which keeps the reduced system symmetric and gives the +// same answer as solving the free block alone. +// +// `values` may be empty, which means every constrained DOF is fixed to zero. +void apply_essential_bc(Sparse& K, std::vector& F, const std::vector& fixed, + const std::vector& values) { const int n = static_cast(F.size()); + // Column pass first: each free row still holds its original entries, so the + // known part of the product can be moved onto F before the columns go. for (int i = 0; i < n; ++i) { + if (fixed[i] != 0) continue; + auto& row = K.rows[i]; + if (!values.empty()) + for (const auto& [j, v] : row) + if (fixed[j] != 0 && values[j] != 0.0) F[i] -= v * values[j]; + row.erase(std::remove_if(row.begin(), row.end(), + [&](const std::pair& e) { + return fixed[e.first] != 0; + }), + row.end()); + } + // Row pass: a constrained DOF becomes the identity equation u_i = g_i. + for (int i = 0; i < n; ++i) if (fixed[i] != 0) { K.rows[i].assign(1, {i, 1.0}); - F[i] = 0.0; - } else { - auto& row = K.rows[i]; - row.erase(std::remove_if(row.begin(), row.end(), - [&](const std::pair& e) { - return fixed[e.first] != 0; - }), - row.end()); + F[i] = values.empty() ? 0.0 : values[i]; } +} + +// Fold `prescribed` into the (fixed, values) pair the elimination takes. A DOF +// already in `fixed` is upgraded to its prescribed value — the same "prescribed +// wins over a plain clamp on the same face" rule the solid path gets from MFEM. +// Two prescribed values for one DOF are a modelling contradiction, not something +// to pick a winner for. +void add_prescribed(const std::vector>& prescribed, int nDof, + const char* what, std::vector& fixed, std::vector& values) { + if (prescribed.empty()) return; + values.assign(nDof, 0.0); + std::vector seen(nDof, 0); + for (const auto& [dof, val] : prescribed) { + if (dof < 0 || dof >= nDof) + throw std::runtime_error(std::string(what) + ": prescribed DOF out of range"); + if (seen[dof] != 0 && values[dof] != val) + throw std::runtime_error( + std::string(what) + ": node " + std::to_string(dof / 6) + " component " + + std::to_string(dof % 6) + " is prescribed twice, to " + + std::to_string(values[dof]) + " and " + std::to_string(val) + + " — a DOF cannot be driven to two displacements at once."); + seen[dof] = 1; + fixed[dof] = 1; + values[dof] = val; } } @@ -611,7 +652,9 @@ ShellResult solve_shell_core(const ShellInput& in) { if (d < 0 || d >= nDof) throw std::runtime_error("shell: fixed DOF out of range"); fixed[d] = 1; } - apply_homogeneous_bc(K, F, fixed); + std::vector values; + add_prescribed(in.prescribed_dofs, nDof, "shell", fixed, values); + apply_essential_bc(K, F, fixed, values); ShellResult res = cg_solve(K, F); if (res.dofs.empty()) res.dofs.assign(nDof, 0.0); @@ -978,6 +1021,7 @@ void resolve_constraint_chains(Rbe3Constraints& C, int nDof) { // coexist in full (WASM heap headroom). ShellResult solve_reduced_system(Sparse& K, const std::vector& F, const std::vector& fixed, + const std::vector& values, const Rbe3Constraints& C, int nDof) { std::vector red(nDof, -1); int nIndep = 0; @@ -1022,27 +1066,38 @@ ShellResult solve_reduced_system(Sparse& K, const std::vector& F, if (F[i] != 0.0) for (const auto& [pi, ci] : expand(i)) Fr[red[pi]] += ci * F[i]; + // An independent DOF maps one-to-one onto its reduced column (coefficient 1), + // so a prescribed value carries over unchanged. A dependent one has no column + // of its own and is refused below, which is why no transform is needed here. std::vector fr(nIndep, 0); + std::vector vr; + if (!values.empty()) vr.assign(nIndep, 0.0); for (int i = 0; i < nDof; ++i) { if (fixed[i] == 0) continue; - // A fixed DOF on a dependent (RBE3 distributing-coupling reference) node - // has no reduced-system column of its own — its motion is a weighted - // average of the coupled solid nodes — so it cannot be constrained here. - // Silently skipping it (the old `&& C.dep[i] == 0` guard) dropped the - // user's constraint and still let CG converge on an under-restrained - // model (issue #377). Refuse loudly instead. - if (C.dep[i] != 0) + // A fixed or prescribed DOF on a dependent (RBE3 distributing-coupling + // reference) node has no reduced-system column of its own — its motion is + // a weighted average of the coupled solid nodes — so it cannot be + // constrained here. Silently skipping it (the old `&& C.dep[i] == 0` + // guard) dropped the user's constraint and still let CG converge on an + // under-restrained model (issue #377). Refuse loudly instead. Driving such + // a node is refused for the same reason it cannot be clamped: the + // coupling, not the BC, already dictates the motion (KOF-210). + if (C.dep[i] != 0) { + const bool driven = !values.empty() && values[i] != 0.0; throw std::runtime_error( - "solve_reduced_system: fixed DOF " + std::to_string(i) + " (node " + - std::to_string(i / 6) + ", component " + std::to_string(i % 6) + + std::string("solve_reduced_system: ") + (driven ? "prescribed" : "fixed") + + " DOF " + std::to_string(i) + " (node " + std::to_string(i / 6) + + ", component " + std::to_string(i % 6) + ") lies on a coupling-dependent node — its motion is governed by " "the RBE3 distributing coupling to the solid, so a direct " "constraint on it cannot be honoured and would otherwise be " "silently dropped. Constrain the coupled solid node(s) instead of " "the shell coupling reference node."); + } fr[red[i]] = 1; + if (!vr.empty()) vr[red[i]] = values[i]; } - apply_homogeneous_bc(Kr, Fr, fr); + apply_essential_bc(Kr, Fr, fr, vr); ShellResult rr = cg_solve(Kr, Fr); if (rr.dofs.empty()) rr.dofs.assign(nIndep, 0.0); @@ -1103,6 +1158,8 @@ ShellResult solve_solid_shell_core(const CoupledInput& in) { if (d < 0 || d >= nDof) throw std::runtime_error("coupled: fixed DOF out of range"); fixed[d] = 1; } + std::vector values; + add_prescribed(in.prescribed_dofs, nDof, "coupled", fixed, values); // A coupling reference node's rotations must NOT be auto-fixed, whichever way // the coupling points. For a distributing coupling all six of its DOFs are // eliminated (they become the RBE3 average of its targets) and a fixed @@ -1128,7 +1185,7 @@ ShellResult solve_solid_shell_core(const CoupledInput& in) { has_rotation[n] = (is_shell[n] != 0 || is_coupling_ref[n] != 0) ? 1 : 0; Rbe3Constraints constraints = build_rbe3_constraints(in, nDof, has_rotation); resolve_constraint_chains(constraints, nDof); - return solve_reduced_system(K, F, fixed, constraints, nDof); + return solve_reduced_system(K, F, fixed, values, constraints, nDof); } // ── Stress recovery ─────────────────────────────────────────────────────────── diff --git a/engine/cpp/shell_core.h b/engine/cpp/shell_core.h index 6d54cec7..0f9ceab8 100644 --- a/engine/cpp/shell_core.h +++ b/engine/cpp/shell_core.h @@ -37,6 +37,10 @@ struct ShellInput { double young = 0.0; // Young's modulus E double poisson = 0.0; // Poisson ratio ν std::vector fixed_dofs; // global DOF indices constrained to zero + // Global DOF index → prescribed value: an INHOMOGENEOUS essential BC (u = g), + // eliminated alongside fixed_dofs rather than instead of it. A DOF listed in + // both takes the prescribed value; two conflicting values for one DOF throw. + std::vector> prescribed_dofs; std::vector> loads; // global DOF index → force/moment }; @@ -148,6 +152,10 @@ struct CoupledInput { std::vector thicknesses; // optional per-triangle thickness std::vector couplings; std::vector fixed_dofs; // global DOF (6·node+comp) fixed to zero + // Global DOF → prescribed value (inhomogeneous essential BC), as ShellInput. + // A prescribed DOF must be INDEPENDENT: like a fixed one, it cannot sit on a + // coupling-dependent node, whose motion the reduction already governs. + std::vector> prescribed_dofs; std::vector> loads; // global DOF → force/moment }; diff --git a/engine/cpp/solve_coupled.cpp b/engine/cpp/solve_coupled.cpp index bc43df10..e53ef5e3 100644 --- a/engine/cpp/solve_coupled.cpp +++ b/engine/cpp/solve_coupled.cpp @@ -187,16 +187,29 @@ val solve_coupled(const val& mesh, const val& coupling, const val& bcs, } in.fixed_dofs = i32_vector(bcs["fixed_dofs"], "bcs.fixed_dofs"); + // Inhomogeneous essential BCs: prescribed_dofs[k] is driven to + // prescribed_vals[k]. Optional — a model with none omits both. + val pdofs_js = bcs["prescribed_dofs"]; + if (!pdofs_js.isUndefined() && !pdofs_js.isNull()) { + std::vector pdofs = i32_vector(pdofs_js, "bcs.prescribed_dofs"); + std::vector pvals = f64_vector(bcs["prescribed_vals"], "bcs.prescribed_vals"); + if (pdofs.size() != pvals.size()) + return error_result("coupled: " + std::to_string(pdofs.size()) + + " bcs.prescribed_dofs but " + std::to_string(pvals.size()) + + " bcs.prescribed_vals — one value per prescribed DOF"); + for (size_t i = 0; i < pdofs.size(); ++i) + in.prescribed_dofs.emplace_back(pdofs[i], pvals[i]); + } std::vector load_dofs = i32_vector(bcs["load_dofs"], "bcs.load_dofs"); std::vector load_vals = f64_vector(bcs["load_vals"], "bcs.load_vals"); for (size_t i = 0; i < load_dofs.size(); ++i) in.loads.emplace_back(load_dofs[i], load_vals[i]); printf("[coupled] solid %zu tets (%zu material%s), shell %zu tris, %zu couplings, " - "%zu fixed, %zu loads\n", + "%zu fixed, %zu prescribed, %zu loads\n", tets.size() / 4, solidE.size(), solidE.size() == 1 ? "" : "s", in.triangles.size() / 3, in.couplings.size(), in.fixed_dofs.size(), - in.loads.size()); + in.prescribed_dofs.size(), in.loads.size()); fflush(stdout); try { diff --git a/engine/cpp/solve_coupled.h b/engine/cpp/solve_coupled.h index a492cfa7..cb36c53a 100644 --- a/engine/cpp/solve_coupled.h +++ b/engine/cpp/solve_coupled.h @@ -16,8 +16,11 @@ // coupling: {ref: Int32Array, offsets: Int32Array, solid: Int32Array} // CSR-style distributing couplings: coupling ref-node k ties to // solid[offsets[k]..offsets[k+1]). -// bcs: {fixed_dofs: Int32Array, load_dofs: Int32Array, load_vals: Float64Array} +// bcs: {fixed_dofs: Int32Array, load_dofs: Int32Array, load_vals: Float64Array, +// prescribed_dofs?: Int32Array, prescribed_vals?: Float64Array} // DOF indices are 6·node+component (0..5 = u,v,w,θx,θy,θz). +// fixed_dofs are constrained to ZERO; prescribed_dofs[k] is driven +// to prescribed_vals[k] (inhomogeneous essential BC). // mat_json: {solid:{young_modulus,poisson_ratio}, shell:{young_modulus,poisson_ratio}} // Returns {displacements: Float64Array} (three translations per node) or {error}. emscripten::val solve_coupled(const emscripten::val& mesh, const emscripten::val& coupling, diff --git a/engine/cpp/solve_shell.cpp b/engine/cpp/solve_shell.cpp index 0e590c4e..caf18116 100644 --- a/engine/cpp/solve_shell.cpp +++ b/engine/cpp/solve_shell.cpp @@ -65,6 +65,26 @@ void add_fixed_vertices(const val& fv_js, int n_nodes, std::vector& fixed) } } +// prescribed_dofs pins one component of a vertex to a NON-ZERO value — an +// inhomogeneous Dirichlet condition (a prescribed-displacement drive, or a +// prescribed rotation on a shell node, which carries all six DOFs). Mirrors +// solve_linear_elastic's `prescribed_dofs` contract, extended from components +// 0..2 to 0..5 because a shell node has rotational DOFs to prescribe. +void add_prescribed_dofs(const val& pdofs_js, int n_nodes, + std::vector>& prescribed) { + if (pdofs_js.isUndefined() || pdofs_js.isNull()) + return; + unsigned n = pdofs_js["length"].as(); + for (unsigned i = 0; i < n; ++i) { + val entry = pdofs_js[i]; + int v = entry["vertex"].as(); + int d = entry["dof"].as(); + kofem::bc::require_valid_vertex(v, n_nodes, "add_prescribed_dofs"); + kofem::bc::require_valid_shell_component(d, v, "add_prescribed_dofs"); + prescribed.emplace_back(6 * v + d, entry["value"].as()); + } +} + // point_loads: force [fx,fy,fz] → DOFs 0..2, optional moment [mx,my,mz] → 3..5. void add_point_loads(const val& loads_js, int n_nodes, std::vector>& loads) { @@ -134,12 +154,13 @@ val solve_shell(const val& mesh, const std::string& mat_json, const std::string& val bcs = parse_json(bcs_json); add_fixed_vertices(bcs["fixed_vertices"], n_nodes, in.fixed_dofs); add_fixed_dofs(bcs["fixed_dofs"], n_nodes, in.fixed_dofs); + add_prescribed_dofs(bcs["prescribed_dofs"], n_nodes, in.prescribed_dofs); add_point_loads(bcs["point_loads"], n_nodes, in.loads); printf("[shell] solve: %d nodes, %d triangles, t=%g, E=%g, nu=%g; " - "%zu fixed DOFs, %zu loads\n", + "%zu fixed DOFs, %zu prescribed DOFs, %zu loads\n", n_nodes, (int)(in.triangles.size() / 3), in.thickness, in.young, - in.poisson, in.fixed_dofs.size(), in.loads.size()); + in.poisson, in.fixed_dofs.size(), in.prescribed_dofs.size(), in.loads.size()); fflush(stdout); kofem::shell::ShellResult r; diff --git a/engine/cpp/solve_shell.h b/engine/cpp/solve_shell.h index faa510c9..96b234ab 100644 --- a/engine/cpp/solve_shell.h +++ b/engine/cpp/solve_shell.h @@ -15,8 +15,11 @@ // object of flat typed arrays {vertices: Float64Array, triangles: Int32Array}. // `mat_json` is {young_modulus, poisson_ratio, thickness}. `bcs_json` is // {fixed_vertices?: int[], fixed_dofs?: [{vertex, dofs:int[]}], +// prescribed_dofs?: [{vertex, dof, value}], // point_loads?: [{vertex, force:[fx,fy,fz], moment?:[mx,my,mz]}]} where DOF -// components are 0..5 = (u,v,w,θx,θy,θz). Returns {displacements: Float64Array} +// components are 0..5 = (u,v,w,θx,θy,θz). prescribed_dofs states an +// INHOMOGENEOUS essential BC (u = value) — fixed_* always mean u = 0. +// Returns {displacements: Float64Array} // (three translations per node, node order) or {error: string} on bad input. emscripten::val solve_shell(const emscripten::val& mesh, const std::string& mat_json, diff --git a/engine/tests/shell_validation.cpp b/engine/tests/shell_validation.cpp index d7b4be35..bfc9849c 100644 --- a/engine/tests/shell_validation.cpp +++ b/engine/tests/shell_validation.cpp @@ -402,6 +402,149 @@ bool coupled_fixed_dependent_throws() { return false; } +// ── Prescribed (inhomogeneous) Dirichlet conditions, KOF-210 ───────────────── +// +// Uniaxial extension of a membrane strip [0,L]x[0,b] in the z=0 plane, driven +// EITHER by a prescribed edge displacement ux = delta on x = L, or by the edge +// force that produces the same strain. Symmetry rollers ux=0 on x=0 and uy=0 on +// y=0; bending DOFs pinned so this is the pure membrane (CST) problem. The exact +// solution is the linear field +// ux = eps*x, uy = -nu*eps*y, eps = delta/L, +// which the CST reproduces exactly, so both drives must land on it and on each +// other. Driving by displacement is the discriminator: an implementation that +// pins a prescribed DOF to zero returns the all-zero field instead. +struct StripResult { + double ux_end; // mean ux on the x = L edge + double uy_side; // uy at the (L, b) corner — the free Poisson contraction +}; + +StripResult strip_uniaxial(bool by_displacement, double delta) { + const double L = 1.0, b = 0.25, t = 0.01, E = 210e9, nu = 0.3; + const int nx = 8, ny = 2; + auto id = [&](int i, int j) { return i * (ny + 1) + j; }; + std::vector V; + std::vector Tr; + for (int i = 0; i <= nx; ++i) + for (int j = 0; j <= ny; ++j) { + V.push_back(L * i / nx); V.push_back(b * j / ny); V.push_back(0.0); + } + for (int i = 0; i < nx; ++i) + for (int j = 0; j < ny; ++j) { + const int a0 = id(i,j), a1 = id(i+1,j), a2 = id(i+1,j+1), a3 = id(i,j+1); + Tr.push_back(a0); Tr.push_back(a1); Tr.push_back(a2); + Tr.push_back(a0); Tr.push_back(a2); Tr.push_back(a3); + } + ShellInput in; + in.vertices = V; in.triangles = Tr; in.thickness = t; in.young = E; in.poisson = nu; + const int nNodes = (nx + 1) * (ny + 1); + for (int nd = 0; nd < nNodes; ++nd) + for (int c : {2, 3, 4, 5}) in.fixed_dofs.push_back(6 * nd + c); // membrane only + for (int j = 0; j <= ny; ++j) in.fixed_dofs.push_back(6 * id(0, j) + 0); // x=0: ux=0 + for (int i = 0; i <= nx; ++i) in.fixed_dofs.push_back(6 * id(i, 0) + 1); // y=0: uy=0 + if (by_displacement) { + for (int j = 0; j <= ny; ++j) in.prescribed_dofs.emplace_back(6 * id(nx, j) + 0, delta); + } else { + // Trapezoidal edge tractions equivalent to sigma = E*delta/L. + const double Ftot = E * (delta / L) * b * t; + for (int j = 0; j <= ny; ++j) { + const double wgt = (j == 0 || j == ny) ? 0.5 : 1.0; + in.loads.emplace_back(6 * id(nx, j) + 0, Ftot * wgt / ny); + } + } + ShellResult r = solve_shell_core(in); + double u = 0.0; + for (int j = 0; j <= ny; ++j) u += r.dofs[6 * static_cast(id(nx, j))]; + return {u / (ny + 1), r.dofs[6 * static_cast(id(nx, ny)) + 1]}; +} + +// Prescribed displacement through the COUPLED assembler and its RBE3 reduction: +// a shell cantilever whose root is distributing-coupled to three anchor nodes, +// with the anchors driven uz = delta and no load at all. The anchors move as a +// body, so the RBE3 average hands the root a pure translation and the whole +// shell must ride along rigidly — tip w = delta. A prescribed value dropped +// anywhere between the input and the reduced system leaves the tip at zero. +double coupled_prescribed_rigid_translation(double delta) { + const double L = 2.0, b = 0.3, t = 0.01, E = 2.1e11, nu = 0.3; + const int nx = 20, ny = 4; + auto id = [&](int i, int j) { return i * (ny + 1) + j; }; + std::vector V; std::vector Tr; std::vector root; + for (int i = 0; i <= nx; ++i) + for (int j = 0; j <= ny; ++j) { V.push_back(L * i / nx); V.push_back(b * j / ny); V.push_back(0); } + for (int i = 0; i < nx; ++i) + for (int j = 0; j < ny; ++j) { + const int a = id(i,j), c = id(i+1,j), d = id(i+1,j+1), e = id(i,j+1); + Tr.push_back(a); Tr.push_back(c); Tr.push_back(d); Tr.push_back(a); Tr.push_back(d); Tr.push_back(e); + } + for (int j = 0; j <= ny; ++j) root.push_back(id(0, j)); + const int aBase = (nx + 1) * (ny + 1); + V.insert(V.end(), {-0.1, 0.0, 0.0, -0.1, b, 0.0, -0.1, b / 2, 0.1}); // 3 anchors + CoupledInput in; + in.n_nodes = aBase + 3; in.vertices = V; in.triangles = Tr; + in.shell_young = E; in.shell_poisson = nu; in.thickness = t; + for (int rn : root) { Coupling cp; cp.ref_node = rn; cp.solid_nodes = {aBase, aBase + 1, aBase + 2}; in.couplings.push_back(cp); } + for (int aa = 0; aa < 3; ++aa) { + in.fixed_dofs.push_back(6 * (aBase + aa) + 0); + in.fixed_dofs.push_back(6 * (aBase + aa) + 1); + in.prescribed_dofs.emplace_back(6 * (aBase + aa) + 2, delta); // the drive + } + ShellResult r = solve_solid_shell_core(in); + double w = 0.0; + for (int j = 0; j <= ny; ++j) w += r.dofs[6 * static_cast(id(nx, j)) + 2]; + return w / (ny + 1); +} + +// A prescribed DOF on a distributing-coupling REFERENCE node is refused for the +// same reason a fixed one is (issue #377): the coupling already governs that +// node's motion, so the drive could only be dropped silently. Same model as +// coupled_fixed_dependent_throws, with the offending constraint driven instead +// of clamped. +bool coupled_prescribed_dependent_throws() { + const double L = 2.0, b = 0.3, t = 0.01, E = 2.1e11, nu = 0.3; + const int nx = 8, ny = 2; + auto id = [&](int i, int j) { return i * (ny + 1) + j; }; + std::vector V; std::vector Tr; std::vector root; + for (int i = 0; i <= nx; ++i) + for (int j = 0; j <= ny; ++j) { V.push_back(L * i / nx); V.push_back(b * j / ny); V.push_back(0); } + for (int i = 0; i < nx; ++i) + for (int j = 0; j < ny; ++j) { + const int a = id(i,j), c = id(i+1,j), d = id(i+1,j+1), e = id(i,j+1); + Tr.push_back(a); Tr.push_back(c); Tr.push_back(d); Tr.push_back(a); Tr.push_back(d); Tr.push_back(e); + } + for (int j = 0; j <= ny; ++j) root.push_back(id(0, j)); + const int aBase = (nx + 1) * (ny + 1); + V.insert(V.end(), {-0.1, 0.0, 0.0, -0.1, b, 0.0, -0.1, b / 2, 0.1}); + CoupledInput in; + in.n_nodes = aBase + 3; in.vertices = V; in.triangles = Tr; + in.shell_young = E; in.shell_poisson = nu; in.thickness = t; + for (int rn : root) { Coupling cp; cp.ref_node = rn; cp.solid_nodes = {aBase, aBase + 1, aBase + 2}; in.couplings.push_back(cp); } + for (int aa = 0; aa < 3; ++aa) for (int c = 0; c < 3; ++c) in.fixed_dofs.push_back(6 * (aBase + aa) + c); + in.prescribed_dofs.emplace_back(6 * root[0] + 0, 1e-3); + try { + solve_solid_shell_core(in); + } catch (const std::exception&) { + return true; + } + return false; +} + +// Two different prescribed values on one DOF is a contradiction, not a +// precedence question — the core must say so rather than pick. +bool shell_conflicting_prescribed_throws() { + std::vector V; std::vector T; + plate_mesh(1.0, 2, V, T); + ShellInput in; + in.vertices = V; in.triangles = T; in.thickness = 0.01; in.young = 1e7; in.poisson = 0.3; + for (int nd = 0; nd < 9; ++nd) for (int c = 0; c < 6; ++c) in.fixed_dofs.push_back(6 * nd + c); + in.prescribed_dofs.emplace_back(6 * 4 + 2, 1e-3); + in.prescribed_dofs.emplace_back(6 * 4 + 2, 2e-3); + try { + solve_shell_core(in); + } catch (const std::exception&) { + return true; + } + return false; +} + } // namespace int main() { @@ -454,6 +597,39 @@ int main() { threw ? "rejected as expected" : "did NOT throw"); } + printf("Prescribed (non-zero) Dirichlet conditions (KOF-210):\n"); + { + const double delta = 1e-4, Lstrip = 1.0, bstrip = 0.25, nu_s = 0.3; + const StripResult pres = strip_uniaxial(true, delta); + const StripResult forced = strip_uniaxial(false, delta); + // THE DISCRIMINATOR: the prescribed edge actually reaches delta. Pinned + // to zero (the old behaviour on every shell path) this reads 0. + check(failures, "shell-prescribed-extension", pres.ux_end, delta, 0.01); + // The Poisson contraction is a SOLVED unknown on the driven edge, so it + // also proves the edge was not over-constrained into a full clamp. + check(failures, "shell-prescribed-poisson", pres.uy_side, + -nu_s * (delta / Lstrip) * bstrip, 0.1); + // Displacement drive and the equivalent force drive are the same + // boundary-value problem — they agree only if the known column + // contribution K[:,p]*g really moved onto the right-hand side. + check(failures, "shell-prescribed-vs-force", pres.ux_end, forced.ux_end, 0.01); + } + { + const double delta = 1e-3; + check(failures, "coupled-prescribed-rigid", + coupled_prescribed_rigid_translation(delta), delta, 0.01); + const bool threw = coupled_prescribed_dependent_throws(); + if (!threw) ++failures; + printf(" [%s] %-28s %s\n", threw ? "PASS" : "FAIL", + "prescribed-on-dependent throws", + threw ? "rejected as expected" : "did NOT throw"); + const bool conflict = shell_conflicting_prescribed_throws(); + if (!conflict) ++failures; + printf(" [%s] %-28s %s\n", conflict ? "PASS" : "FAIL", + "conflicting prescribed throws", + conflict ? "rejected as expected" : "did NOT throw"); + } + printf(failures != 0 ? "\n%d check(s) FAILED\n" : "\nall checks passed\n", failures); return failures != 0 ? 1 : 0; } diff --git a/examples/validation/README.md b/examples/validation/README.md index ee1b03cf..bf91ef38 100644 --- a/examples/validation/README.md +++ b/examples/validation/README.md @@ -82,3 +82,26 @@ node examples/validation/prescribed-displacement.test.mjs Like the single-DOF test, this needs the `prescribed_dofs` support compiled into `engine/cpp/engine.cpp`, so it only passes against a **freshly built** WASM and fails on purpose on the committed binary until `scripts/build-wasm.sh` runs. + +## `shell-prescribed-displacement.test.mjs` — non-zero Dirichlet on shells + +The shell counterpart of the test above, guarding KOF-210. The solid path got +inhomogeneous Dirichlet conditions from MFEM; the shell formulation's boundary +conditions were homogeneous-only, so a prescribed displacement on a shell model +either refused to run or — on the auto-shell and coupled paths — silently pinned +the value to zero and returned an all-zero field that looks converged. Thin +walls become shells automatically at mesh time, so this was easy to hit without +choosing a shell idealisation. + +A flat strip is stretched in its own plane by a prescribed edge displacement, +with **no applied load**. The driven edge must reach ux ≈ δ, midspan must land +on the linear field δ/2, and the driven edge must still contract under Poisson's +effect — proving it was driven, not clamped. + +```bash +node examples/validation/shell-prescribed-displacement.test.mjs +``` + +The element math behind it is also covered natively, without a WASM build, by +`bash scripts/test-shell.sh` (the `shell-prescribed-*` and `coupled-prescribed-*` +checks in `engine/tests/shell_validation.cpp`). diff --git a/examples/validation/lib/solver.mjs b/examples/validation/lib/solver.mjs index 95c18e14..94456423 100644 --- a/examples/validation/lib/solver.mjs +++ b/examples/validation/lib/solver.mjs @@ -35,15 +35,7 @@ const pkg = join(here, "../../../web/src/wasm/pkg"); * Returns { displacements:number[] (3/node), von_mises:number[] (1/elem) }. */ export async function loadSolver() { - const wasmBinary = readFileSync(join(pkg, "kofem_wasm_emcc.wasm")).buffer; - const { default: createModule } = await import( - join(pkg, "kofem_wasm_emcc.js") - ); - const Module = await createModule({ - wasmBinary, - print: () => {}, - printErr: () => {}, - }); + const Module = await loadModule(); return function solve(mesh, material, bcs, order = 1) { // The engine takes the mesh as flat typed arrays (issue #166); the @@ -75,3 +67,52 @@ export async function loadSolver() { }; }; } + +// The raw Emscripten module, shared by the solve wrappers below. +async function loadModule() { + const wasmBinary = readFileSync(join(pkg, "kofem_wasm_emcc.wasm")).buffer; + const { default: createModule } = await import( + join(pkg, "kofem_wasm_emcc.js") + ); + return createModule({ wasmBinary, print: () => {}, printErr: () => {} }); +} + +/** + * Initialise the engine and return a solveShell() closure over the Kirchhoff + * flat-facet shell entry point (solve_shell). Same nested-tuple → typed-array + * flattening as loadSolver, for the shell's triangle surface mesh: + * mesh: { vertices:[[x,y,z]...], triangles:[[a,b,c]...], + * thicknesses?:[t...] } (per-facet thickness) + * material: { young_modulus, poisson_ratio, thickness } + * bcs: { fixed_vertices:[v...], + * fixed_dofs:[{vertex, dofs:[0..5,...]}...], // u = 0 + * prescribed_dofs:[{vertex, dof:0..5, value}...], // u = value + * point_loads:[{vertex, force:[fx,fy,fz], moment?:[mx,my,mz]}...] } + * Shell DOF components are 0..5 = (u,v,w,θx,θy,θz). + * Returns { displacements:number[] (3/node), von_mises:number[] (1/triangle) }. + */ +export async function loadShellSolver() { + const Module = await loadModule(); + + return function solveShell(mesh, material, bcs) { + const result = Module.solve_shell( + { + vertices: Float64Array.from(mesh.vertices.flat()), + triangles: Int32Array.from(mesh.triangles.flat()), + thicknesses: Float64Array.from(mesh.thicknesses ?? []), + }, + JSON.stringify(material), + JSON.stringify({ + fixed_vertices: bcs.fixed_vertices ?? [], + fixed_dofs: bcs.fixed_dofs ?? [], + prescribed_dofs: bcs.prescribed_dofs ?? [], + point_loads: bcs.point_loads ?? [], + }), + ); + if ("error" in result) throw new Error(result.error); + return { + displacements: Array.from(result.displacements), + von_mises: Array.from(result.von_mises), + }; + }; +} diff --git a/examples/validation/shell-prescribed-displacement.test.mjs b/examples/validation/shell-prescribed-displacement.test.mjs new file mode 100644 index 00000000..94bf713f --- /dev/null +++ b/examples/validation/shell-prescribed-displacement.test.mjs @@ -0,0 +1,168 @@ +// SPDX-FileCopyrightText: 2026 Michael Kofler +// SPDX-License-Identifier: AGPL-3.0-or-later + +// Regression test for non-zero prescribed displacements on SHELL models (KOF-210). +// Shell counterpart of prescribed-displacement.test.mjs, which covers the solid +// (MFEM) path. +// +// THE BUG THIS GUARDS: the shell formulation's boundary conditions used to be +// homogeneous-only — ShellInput.fixed_dofs are "constrained to zero" and there +// was no inhomogeneous counterpart to the solid path's prescribed_dofs. A model +// driven by a prescribed displacement therefore either refused to run (explicit +// CTRIA3) or, worse, silently pinned the value to zero and returned an all-zero +// field that looks like a converged answer (auto-shell / coupled). Thin-walled +// parts land on the shell path automatically, so this was easy to hit without +// ever choosing a shell idealisation. +// +// Setup: in-plane (membrane) extension of a flat strip [0,L]x[0,b] in the z = 0 +// plane, driven by displacement only, with no applied load. Symmetry rollers +// Ux = 0 on x = 0 and Uy = 0 on y = 0, a NON-ZERO prescribed Ux = δ on the x = L +// edge, and the out-of-plane DOFs (w, θx, θy) plus the stiffness-free drilling +// DOF θz pinned everywhere so this is the pure CST membrane problem. The exact +// solution is the linear field +// ux = εx·x, uy = −ν·εx·y, εx = δ/L, +// independent of E and t (a pure-displacement boundary-value problem), which the +// constant-strain triangle reproduces exactly. The x = L edge must reach ux ≈ δ; +// a binary without shell prescribed_dofs support leaves it at zero. +// +// Runs against the WASM engine the other validation tests use. On a binary built +// before the shell prescribed_dofs support it is EXPECTED to fail — either with +// the old refusal or with an all-zero field. That failure is the signal to +// rebuild with scripts/build-wasm.sh. + +import { loadShellSolver } from "./lib/solver.mjs"; +import { nodesWhere } from "./lib/mesh.mjs"; + +const nu = 0.3; +const E = 210e9; // result is E-independent; any positive value works +const t = 0.01; +const L = 1.0, + b = 0.25; +const delta = 1e-4; // prescribed Ux on the x = L edge +const epsX = delta / L; +const expUySide = -nu * epsX * b; // Poisson contraction at y = b + +// Structured triangle mesh of the strip, two triangles per cell. +const nx = 8, + ny = 2; +const id = (i, j) => i * (ny + 1) + j; +const vertices = []; +for (let i = 0; i <= nx; i++) + for (let j = 0; j <= ny; j++) vertices.push([(L * i) / nx, (b * j) / ny, 0]); +const triangles = []; +for (let i = 0; i < nx; i++) + for (let j = 0; j < ny; j++) { + triangles.push([id(i, j), id(i + 1, j), id(i + 1, j + 1)]); + triangles.push([id(i, j), id(i + 1, j + 1), id(i, j + 1)]); + } + +// Zero-Dirichlet set, unioned per vertex: the out-of-plane and drilling DOFs +// everywhere (a flat plate loaded in its own plane has no bending action, and +// DKT+CST carries no drilling stiffness, so those would otherwise be singular), +// plus the two symmetry rollers. +const zeroByVertex = new Map(); +const addZero = (ids, dofs) => + ids.forEach((i) => { + if (!zeroByVertex.has(i)) zeroByVertex.set(i, new Set()); + for (const d of dofs) zeroByVertex.get(i).add(d); + }); +addZero( + vertices.map((_v, i) => i), + [2, 3, 4, 5], +); +addZero( + nodesWhere(vertices, (x) => x <= 1e-9), + [0], +); // x = 0 → Ux = 0 +addZero( + nodesWhere(vertices, (x, y) => y <= 1e-9), + [1], +); // y = 0 → Uy = 0 +const fixed_dofs = [...zeroByVertex].map(([vertex, s]) => ({ + vertex, + dofs: [...s].sort((p, q) => p - q), +})); + +// The driving condition: NON-ZERO prescribed Ux on the x = L edge, no load. +const driven = nodesWhere(vertices, (x) => x >= L - 1e-9); +const prescribed_dofs = driven.map((vertex) => ({ + vertex, + dof: 0, + value: delta, +})); + +const solveShell = await loadShellSolver(); +const r = solveShell( + { vertices, triangles }, + { young_modulus: E, poisson_ratio: nu, thickness: t }, + { fixed_vertices: [], fixed_dofs, prescribed_dofs, point_loads: [] }, +); + +const d = (v, c) => r.displacements[v * 3 + c]; +const uxEnd = driven.reduce((s, v) => s + d(v, 0), 0) / driven.length; // ≈ δ +// Discriminator vertex: on the driven edge AND at y = b, so it is free to +// contract under Poisson — the (L, b) corner. +const corner = nodesWhere( + vertices, + (x, y) => x >= L - 1e-9 && y >= b - 1e-9, +)[0]; +// Midspan vertex on the symmetry line y = 0, where the linear field gives δ/2. +const mid = nodesWhere( + vertices, + (x, y) => Math.abs(x - L / 2) <= 1e-9 && y <= 1e-9, +)[0]; + +const checks = []; +const check = (name, ok, detail) => checks.push({ name, ok, detail }); + +const finite = r.displacements.every(Number.isFinite); +check("shell solve produced finite displacements", finite, ""); + +if (finite) { + // THE DISCRIMINATOR: the prescribed value actually reaches the driven edge. + // A binary without shell prescribed_dofs pins it to zero, so "within 1% of δ" + // cleanly separates applied from discarded. + check( + "prescribed Ux is applied on the driven edge (ux ≈ δ, NOT zero)", + Math.abs((uxEnd - delta) / delta) < 0.01, + `ux=${uxEnd.toExponential(3)} vs δ=${delta.toExponential(3)}`, + ); + + // The CST reproduces the linear extension field exactly, so midspan lands on + // δ/2 — proof the whole field is driven, not just the constrained edge. + check( + "midspan follows the linear extension field (ux ≈ δ/2)", + Math.abs((d(mid, 0) - delta / 2) / (delta / 2)) < 0.01, + `ux=${d(mid, 0).toExponential(3)} vs ${(delta / 2).toExponential(3)}`, + ); + + // The transverse contraction is a SOLVED unknown, so it proves the driven edge + // is FREE in y rather than clamped: an implementation that folded the + // prescribed BC into a full fix would return uy = 0 here. + check( + "driven edge is FREE in Uy and contracts under Poisson (not over-pinned)", + Math.abs((d(corner, 1) - expUySide) / expUySide) < 0.05, + `uy=${d(corner, 1).toExponential(3)} vs ${expUySide.toExponential(3)}`, + ); +} + +console.log("\nShell prescribed-displacement test (KOF-210)\n"); +let failed = 0; +for (const c of checks) { + console.log( + ` ${c.ok ? "PASS" : "FAIL"} ${c.name}${c.detail ? " — " + c.detail : ""}`, + ); + if (!c.ok) failed++; +} + +if (failed) { + console.error( + `\n${failed} check(s) FAILED. If this is a stale WASM binary, rebuild it with` + + `\n scripts/build-wasm.sh` + + `\nso shell_core's prescribed_dofs (inhomogeneous Dirichlet) support is compiled in.\n`, + ); + process.exit(1); +} +console.log( + "\nPASS — non-zero prescribed displacements are honored on shell models.\n", +); diff --git a/web/package.json b/web/package.json index e5c4200e..82cf79f2 100644 --- a/web/package.json +++ b/web/package.json @@ -23,7 +23,7 @@ "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", "pretest": "bun run wasm:fetch", - "test": "bun ../examples/validation/run.mjs && bun ../examples/validation/dof-constraint.test.mjs && bun ../examples/validation/prescribed-displacement.test.mjs && bun ../examples/validation/multiple-loads.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_multi_load.mjs && bun tests/test_shellize_mpc.mjs && bun tests/test_mesh_sizing.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 ../examples/validation/shell-prescribed-displacement.test.mjs && bun ../examples/validation/multiple-loads.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_multi_load.mjs && bun tests/test_shellize_mpc.mjs && bun tests/test_mesh_sizing.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/wasm/pkg/kofem_wasm.d.ts b/web/src/wasm/pkg/kofem_wasm.d.ts index bdf6f7cf..7db6e5f5 100644 --- a/web/src/wasm/pkg/kofem_wasm.d.ts +++ b/web/src/wasm/pkg/kofem_wasm.d.ts @@ -104,9 +104,11 @@ export interface KofemModule { ): StaticSolveResult /** Kirchhoff flat-facet shell solve on a triangle surface mesh. `mat_json` is * `{ young_modulus, poisson_ratio, thickness }`; `bcs_json` is - * `{ fixed_vertices?, fixed_dofs?, point_loads? }` with DOF components - * 0..5 = (u,v,w,θx,θy,θz). Returns three translations per node plus one von - * Mises surface stress per triangle. */ + * `{ fixed_vertices?, fixed_dofs?, prescribed_dofs?, point_loads? }` with DOF + * components 0..5 = (u,v,w,θx,θy,θz). `fixed_*` always mean u = 0; + * `prescribed_dofs: [{vertex, dof, value}]` drives a DOF to a non-zero value + * (inhomogeneous essential BC). Returns three translations per node plus one + * von Mises surface stress per triangle. */ solve_shell( mesh: ShellMesh, mat_json: string, @@ -128,7 +130,10 @@ export interface KofemModule { * 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) + * bcs: { fixed_dofs, load_dofs, load_vals, prescribed_dofs?, + * prescribed_vals? } (DOF = 6·node+component). fixed_dofs are + * constrained to ZERO; prescribed_dofs[k] is driven to + * prescribed_vals[k] (inhomogeneous essential BC). * mat_json: { solid, shell:{young_modulus,poisson_ratio} } — `solid` is * either one material object (every tet uses it) or an ARRAY of * materials selected per tet by `mesh.attributes`. The shell side @@ -149,7 +154,13 @@ export interface KofemModule { dof_mask?: Int32Array relaxation?: number }, - bcs: { fixed_dofs: Int32Array; load_dofs: Int32Array; load_vals: Float64Array }, + bcs: { + fixed_dofs: Int32Array + load_dofs: Int32Array + load_vals: Float64Array + prescribed_dofs?: Int32Array + prescribed_vals?: Float64Array + }, mat_json: string, ): | { diff --git a/web/src/workers/solver.worker.ts b/web/src/workers/solver.worker.ts index 9c0149a5..35dfaa45 100644 --- a/web/src/workers/solver.worker.ts +++ b/web/src/workers/solver.worker.ts @@ -852,22 +852,27 @@ function resolveShellSection( // Essential BCs for the shell solve. Shell nodes carry six DOFs // (u,v,w,θx,θy,θz), so rotational constraints are honoured — unlike the solid // path, which drops them as stiffness-free. A node with all six DOFs fixed -// becomes a fixed_vertices entry, anything partial a fixed_dofs entry. The -// shell solver has no inhomogeneous essential BCs, so a non-zero prescribed -// displacement is a loud error, not a silent pin-to-zero. +// becomes a fixed_vertices entry, anything partial a fixed_dofs entry. +// +// A non-zero prescribed displacement (or rotation) is split off into +// prescribed_dofs exactly as groupDirichlet does for the solid path: the +// fixed_* sets always mean u = 0, so folding a driven DOF into them would pin +// away the value the user asked for (KOF-210). function shellDirichlet(constraints: Constraint[], vid: VertexIndexer) { const dofsByVertex = new Map>(); + const prescribed_dofs: { vertex: number; dof: number; value: number }[] = []; for (const c of constraints) { if (c.dof < 0 || c.dof > 5) throw new Error( `shell solve: constraint on node ${c.nodeId} names DOF ${c.dof} — valid shell DOFs are 0..5`, ); - // eslint-disable-next-line kofem/no-silent-fallback -- a constraint without prescribedValue is a homogeneous fixed BC, i.e. u = 0 by definition - if ((c.prescribedValue ?? 0) !== 0) - throw new Error( - "shell solve: prescribed (non-zero) displacements are not supported for shell models yet", - ); const vertex = vid(c.nodeId, "constraint"); + // eslint-disable-next-line kofem/no-silent-fallback -- a constraint without prescribedValue is a homogeneous fixed BC, i.e. u = 0 by definition + const value = c.prescribedValue ?? 0; + if (value !== 0) { + prescribed_dofs.push({ vertex, dof: c.dof, value }); + continue; + } let dofs = dofsByVertex.get(vertex); if (!dofs) { dofs = new Set(); @@ -881,7 +886,7 @@ function shellDirichlet(constraints: Constraint[], vid: VertexIndexer) { if (dofSet.size === 6) fixed_vertices.push(vertex); else fixed_dofs.push({ vertex, dofs: [...dofSet].sort((a, b) => a - b) }); } - return { fixed_vertices, fixed_dofs }; + return { fixed_vertices, fixed_dofs, prescribed_dofs }; } type ShellNodalLoad = { @@ -1069,7 +1074,10 @@ function handleShellSolve(id: number, payload: SolvePayload) { materials, properties, ); - const { fixed_vertices, fixed_dofs } = shellDirichlet(constraints, vid); + const { fixed_vertices, fixed_dofs, prescribed_dofs } = shellDirichlet( + constraints, + vid, + ); const posOf = (nodeId: number): [number, number, number] => { const vi = vid(nodeId, "surface load face"); return [vertices[3 * vi], vertices[3 * vi + 1], vertices[3 * vi + 2]]; @@ -1083,7 +1091,12 @@ function handleShellSolve(id: number, payload: SolvePayload) { const result = m().solve_shell( { vertices, triangles, thicknesses }, JSON.stringify({ young_modulus: young, poisson_ratio: poisson }), - JSON.stringify({ fixed_vertices, fixed_dofs, point_loads }), + JSON.stringify({ + fixed_vertices, + fixed_dofs, + prescribed_dofs, + point_loads, + }), ); if ("error" in result) throw new Error(result.error); @@ -1213,29 +1226,55 @@ function coupledMaterials( // 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. +// +// A non-zero prescribed value leaves through `prescribed_dofs`/`prescribed_vals` +// (the coupled engine's inhomogeneous essential BCs), never through fixed_dofs, +// which always mean u = 0 (KOF-210). A driven node is NOT given the +// all-translations-clamped rotational treatment: a face pulled to a displacement +// is not thereby built in, and clamping its rotations would over-stiffen it. function coupledFixedDofs( constraints: Constraint[], poolOf: (nodeId: number) => number, isShell: (poolIndex: number) => boolean, isRefPoint: (poolIndex: number) => boolean, -): number[] { +): { + fixed_dofs: number[]; + prescribed_dofs: number[]; + prescribed_vals: number[]; +} { const fixedByPool = new Map>(); + const drivenPools = new Set(); + const prescribed = new Map(); for (const c of constraints) { - // The coupled assembler's essential BCs are homogeneous — shell_core's - // ShellInput.fixed_dofs are "constrained to zero", with no inhomogeneous - // counterpart to the solid path's prescribed_dofs. A non-zero prescribed - // displacement used to be pinned to zero here without a word, which turns a - // displacement-driven model into an unloaded one and returns an all-zero - // field that looks like a converged answer (KOF-216). - // eslint-disable-next-line kofem/no-silent-fallback -- a constraint without prescribedValue is a homogeneous fixed BC, i.e. u = 0 by definition - if ((c.prescribedValue ?? 0) !== 0) - throw new Error( - "Prescribed (non-zero) displacements are not supported on the coupled shell/solid " + - `path yet: node ${c.nodeId} prescribes ${c.prescribedValue} on DOF ${c.dof}. ` + - "Drive this model with a load instead, or solve it without the shell idealisation.", - ); const pi = poolOf(c.nodeId); - if (c.dof > 2 && !isRefPoint(pi)) continue; + // eslint-disable-next-line kofem/no-silent-fallback -- a constraint without prescribedValue is a homogeneous fixed BC, i.e. u = 0 by definition + const value = c.prescribedValue ?? 0; + if (c.dof > 2 && !isRefPoint(pi)) { + // A rotational CLAMP is dropped here by design (see above), but a driven + // rotation has no such substitute — dropping it is the silent pin-to-zero + // this path was fixed for (KOF-210). Say so instead. + if (value !== 0) + throw new Error( + `Node ${c.nodeId} prescribes a rotation (DOF ${c.dof} = ${value}), which this ` + + "model has nowhere to apply: only a coupling reference point carries a " + + "drivable rotational DOF here. Prescribe the translations instead, or " + + "declare the node as a coupling reference point.", + ); + continue; + } + if (value !== 0) { + const dof = 6 * pi + c.dof; + const seen = prescribed.get(dof); + if (seen !== undefined && seen !== value) + throw new Error( + `Node ${c.nodeId} (DOF ${c.dof}) is prescribed to two different displacements, ` + + `${seen} and ${value} — a DOF cannot be driven to both. Remove or merge one ` + + "of the boundary conditions.", + ); + prescribed.set(dof, value); + drivenPools.add(pi); + continue; + } let dofs = fixedByPool.get(pi); if (!dofs) { dofs = new Set(); @@ -1245,7 +1284,12 @@ function coupledFixedDofs( } const fixed_dofs: number[] = []; for (const [pi, dofs] of fixedByPool) { - for (const d of dofs) fixed_dofs.push(6 * pi + d); + // A DOF driven to a value wins over a plain clamp on the same DOF, the same + // way the solid path's prescribed_dofs override an overlapping fixed group. + for (const d of dofs) { + if (prescribed.has(6 * pi + d)) continue; + fixed_dofs.push(6 * pi + d); + } // 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 @@ -1254,13 +1298,20 @@ function coupledFixedDofs( if ( isShell(pi) && !isRefPoint(pi) && + !drivenPools.has(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; + const prescribed_dofs: number[] = []; + const prescribed_vals: number[] = []; + for (const [dof, value] of prescribed) { + prescribed_dofs.push(dof); + prescribed_vals.push(value); + } + return { fixed_dofs, prescribed_dofs, prescribed_vals }; } // Point + surface loads → equivalent nodal forces on the pool. @@ -1485,20 +1536,24 @@ function tryCoupledSolve( const refPoolIndices = new Set(model.refPool.values()); const isRefPoint = (pi: number) => refPoolIndices.has(pi); - const fixed_dofs = coupledFixedDofs( + const { fixed_dofs, prescribed_dofs, prescribed_vals } = 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. + // A clamped or DRIVEN shell rim can sit next to the retained base solid; the + // proximity detector would otherwise couple the very nodes the user constrained + // (the engine refuses a constrained coupling-dependent node, #377/KOF-210). The + // BC wins, for a prescribed displacement exactly as for a clamp. // 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), + dropCouplingsOnFixedNodes(model.coupling, [ + ...fixed_dofs, + ...prescribed_dofs, + ]), buildReferenceCouplings(couplings, (nodeId) => poolOf(nodeId)), ); const { load_dofs, load_vals } = coupledLoads( @@ -1545,6 +1600,8 @@ function tryCoupledSolve( }, { fixed_dofs: Int32Array.from(fixed_dofs), + prescribed_dofs: Int32Array.from(prescribed_dofs), + prescribed_vals: Float64Array.from(prescribed_vals), load_dofs: Int32Array.from(load_dofs), load_vals: Float64Array.from(load_vals), }, @@ -1720,14 +1777,14 @@ function tryPureShellSolve( // 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( + const flat = coupledFixedDofs( constraints, shellOf, () => true, noReferencePoints, ); const dofsByVertex = new Map>(); - for (const d of flatFixed) + for (const d of flat.fixed_dofs) getOrInitDofs(dofsByVertex, Math.floor(d / 6)).add(d % 6); const fixed_vertices: number[] = []; const fixed_dofs: { vertex: number; dofs: number[] }[] = []; @@ -1735,6 +1792,13 @@ function tryPureShellSolve( if (dofSet.size === 6) fixed_vertices.push(vertex); else fixed_dofs.push({ vertex, dofs: [...dofSet].sort((a, b) => a - b) }); } + // Driven DOFs stay separate: solve_shell's fixed_* sets mean u = 0, so a + // prescribed value only survives as an inhomogeneous essential BC (KOF-210). + const prescribed_dofs = flat.prescribed_dofs.map((d, k) => ({ + vertex: Math.floor(d / 6), + dof: d % 6, + value: flat.prescribed_vals[k], + })); const { load_dofs, load_vals } = coupledLoads( loads, @@ -1770,7 +1834,12 @@ function tryPureShellSolve( thicknesses: Float64Array.from(shells.shellThk), }, JSON.stringify({ young_modulus: mat.young, poisson_ratio: mat.poisson }), - JSON.stringify({ fixed_vertices, fixed_dofs, point_loads }), + JSON.stringify({ + fixed_vertices, + fixed_dofs, + prescribed_dofs, + point_loads, + }), ); if ("error" in result) throw new Error(result.error); @@ -2047,17 +2116,21 @@ function handleMixedSolve(id: number, payload: SolvePayload) { }; const refPoolIndices = new Set([...refIds].map((nodeId) => poolOf(nodeId))); const isRefPoint = (pi: number) => refPoolIndices.has(pi); - const fixed_dofs = coupledFixedDofs( + const { fixed_dofs, prescribed_dofs, prescribed_vals } = 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. + // A clamped or driven shell node that also sits within coupling range of the + // solid would be both constrained and a distributing-coupling dependent — the + // engine refuses that (#377/KOF-210). The BC wins; drop the coupling on those + // nodes. const coupling = concatCouplings( - dropCouplingsOnFixedNodes(model.coupling, fixed_dofs), + dropCouplingsOnFixedNodes(model.coupling, [ + ...fixed_dofs, + ...prescribed_dofs, + ]), buildReferenceCouplings(couplings, (nodeId, context) => { const pi = model.poolOfVertex.get(vid(nodeId, context)); if (pi === undefined) @@ -2104,6 +2177,8 @@ function handleMixedSolve(id: number, payload: SolvePayload) { }, { fixed_dofs: Int32Array.from(fixed_dofs), + prescribed_dofs: Int32Array.from(prescribed_dofs), + prescribed_vals: Float64Array.from(prescribed_vals), load_dofs: Int32Array.from(load_dofs), load_vals: Float64Array.from(load_vals), },