Support prescribed (non-zero) displacements on shell models

engine/cpp/shell_core.cpp

@@ -555,20 +555,61 @@ std::array<std::array<double, 3>, 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<double>& F, const std::vector<char>& 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<double>& F, const std::vector<char>& fixed,
+                        const std::vector<double>& values) {
     const int n = static_cast<int>(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<int, double>& 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<int, double>& 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<std::pair<int, double>>& prescribed, int nDof,
+                    const char* what, std::vector<char>& fixed, std::vector<double>& values) {
+    if (prescribed.empty()) return;
+    values.assign(nDof, 0.0);
+    std::vector<char> 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<double> 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<double>& F,
                                  const std::vector<char>& fixed,
+                                 const std::vector<double>& values,
                                  const Rbe3Constraints& C, int nDof) {
     std::vector<int> red(nDof, -1);
     int nIndep = 0;
@@ -1022,27 +1066,38 @@ ShellResult solve_reduced_system(Sparse& K, const std::vector<double>& 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<char> fr(nIndep, 0);
+    std::vector<double> 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<double> 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 ───────────────────────────────────────────────────────────

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<int> 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<std::pair<int, double>> prescribed_dofs;
     std::vector<std::pair<int, double>> loads;  // global DOF index → force/moment
 };
 
@@ -148,6 +152,10 @@ struct CoupledInput {
     std::vector<double> thicknesses;     // optional per-triangle thickness
     std::vector<Coupling> couplings;
     std::vector<int> 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<std::pair<int, double>> prescribed_dofs;
     std::vector<std::pair<int, double>> loads;  // global DOF → force/moment
 };
 

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<int> pdofs = i32_vector(pdofs_js, "bcs.prescribed_dofs");
+        std::vector<double> 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<int> load_dofs = i32_vector(bcs["load_dofs"], "bcs.load_dofs");
     std::vector<double> 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 {

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,

engine/cpp/solve_shell.cpp

@@ -65,6 +65,26 @@ void add_fixed_vertices(const val& fv_js, int n_nodes, std::vector<int>& 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<std::pair<int, double>>& prescribed) {
+    if (pdofs_js.isUndefined() || pdofs_js.isNull())
+        return;
+    unsigned n = pdofs_js["length"].as<unsigned>();
+    for (unsigned i = 0; i < n; ++i) {
+        val entry = pdofs_js[i];
+        int v = entry["vertex"].as<int>();
+        int d = entry["dof"].as<int>();
+        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<double>());
+    }
+}
+
 // 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<std::pair<int, double>>& 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;

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,