Use a 2×2 Gauss rule for linear quadrilaterals (#105)

Linear elements used a 1-point Gauss rule per direction, which
under-integrates the stiffness matrix of bilinear quadrilaterals (and the
convection terms on their sides), so results converged slowly; on the
NAFEMS T4 benchmark a 6 × 10 linear mesh gave 28.27 °C instead of
18.25 °C. Linear elements now use a 2-point rule per direction (exact up
to degree 3), giving 17.95 °C on the same mesh.

Tests: the quadrature unit test checks exactness up to degree 3 for
linear quadrilaterals, and the regression values of the linear cases
(advection-diffusion 1D, heat conduction 2D fin, nonlinear
reaction-diffusion 1D, solidification front 2D) are updated, with the
reason in each README.

Closes #105

Add NAFEMS T4, Couette flow and cantilever beam examples

Node.js examples based on the verification tests, each printing the
result next to its reference value:
- heatConduction2DNafemsT4: NAFEMS T4 benchmark (heat conduction with
  convection), temperature at point E against the 18.3 °C target
  (18.25 °C converged)
- couetteFlow2DCreepingFlow: Couette flow, against the exact linear
  velocity profile with zero pressure
- cantileverBeam1D: cantilever with a tip load, against the closed-form
  tip deflection and rotation

Author: nikoscham

dist/feascript-worker.esm.js

@@ -19,6 +19,13 @@ This example solves a 2D lid-driven cavity flow using the creeping flow solver w
 [tutorial](https://feascript.com/tutorials/lid-driven-cavity-2d-creeping-flow.html) on the
 FEAScript website.
 
+### 2. Couette Flow in a 2D Channel (`couetteFlow2DCreepingFlow.js`)
+
+A channel whose top wall moves with velocity U, with a fixed bottom wall and stress-free (open)
+ends. The exact solution is the linear velocity profile u = U y / H with v = 0 and zero pressure,
+which the Taylor-Hood elements reproduce exactly. The example prints the velocity along the left
+end and the largest deviation from the exact solution.
+
 ## Running the Node.js examples
 
 ### 1. Create `package.json` with ES module support
@@ -39,4 +46,6 @@ npm install feascript
 node lidDrivenCavity2DCreepingFlow.js
 ```
 
+(or `node couetteFlow2DCreepingFlow/couetteFlow2DCreepingFlow.js`)
+
 (for the "Creeping Flow on a Lid-Driven 2D Cavity" example)

dist/feascript-worker.esm.js.map

@@ -0,0 +1,65 @@
+/**
+ * ════════════════════════════════════════════════════════════════
+ *  FEAScript Core Library
+ *  Finite Element Simulation in JavaScript
+ *  Version: 0.3.0 | https://feascript.com
+ *  MIT License © 2023–2026 FEAScript
+ * ════════════════════════════════════════════════════════════════
+ */
+
+// Couette flow: a 4 × 2 channel whose top wall moves with velocity U, with stress-free (open) ends.
+// The exact solution is the linear velocity profile u = U y / H, v = 0, with zero pressure, which
+// Taylor-Hood elements reproduce exactly
+
+// Import FEAScript library
+import { FEAScriptModel, printVersion } from "feascript";
+
+console.log("FEAScript Version:", printVersion);
+
+const wallVelocity = 1; // U
+const channelHeight = 2; // H
+
+// 1. Create the model and select the physics
+const model = new FEAScriptModel();
+model.setModelConfig("creepingFlowScript", {
+  coefficientFunctions: {
+    viscosity: 1, // μ
+  },
+});
+
+// 2. Define the mesh (creeping flow needs quadratic elements)
+model.setMeshConfig({
+  meshDimension: "2D",
+  elementOrder: "quadratic",
+  numElementsX: 4,
+  numElementsY: 2,
+  maxX: 4,
+  maxY: channelHeight,
+});
+
+// 3. Apply the boundary conditions
+model.addBoundaryCondition("0", ["constantVelocity", 0, 0]); // Bottom: fixed wall
+model.addBoundaryCondition("1", ["stressFree"]); // Left: open end
+model.addBoundaryCondition("2", ["constantVelocity", wallVelocity, 0]); // Top: moving wall
+model.addBoundaryCondition("3", ["stressFree"]); // Right: open end
+
+// 4. Solve
+model.setSolverMethod("lusolve");
+const { nodesCoordinates, fields } = model.solve();
+
+// 5. Print the results
+const { nodesYCoordinates } = nodesCoordinates;
+const largestDeviation = Math.max(
+  ...nodesYCoordinates.map((y, i) => Math.abs(fields.u[i] - (wallVelocity * y) / channelHeight)),
+);
+console.log("Velocity u along the left end (x = 0):");
+nodesYCoordinates.forEach((y, i) => {
+  if (nodesCoordinates.nodesXCoordinates[i] === 0) {
+    console.log(
+      `  y = ${y.toFixed(2)}: u = ${fields.u[i].toFixed(4)} (exact ${(wallVelocity * y) / channelHeight})`,
+    );
+  }
+});
+console.log(`Largest deviation from u = U y / H: ${largestDeviation.toExponential(2)}`);
+console.log(`Largest |v|: ${Math.max(...fields.v.map(Math.abs)).toExponential(2)} (exact 0)`);
+console.log(`Largest |p|: ${Math.max(...fields.p.map(Math.abs)).toExponential(2)} (exact 0)`);