feature 1D-PDE-Solver for FEAScript-core (#42)

This commit merges the PR introducing a general-purpose solver for 1D linear PDEs, supporting equations like d/dx(A du/dx) + B du/dx + C u = D with user-defined coefficients and boundary conditions. The implementation is standalone and functional as a prototype. Further integration is needed to align with FEAScript's existing utilities.

Author: Stu-ops

examples/generalFormPDEScript/exampleGeneralFormPDE.js

@@ -0,0 +1,33 @@
+//   ______ ______           _____           _       _     //
+//  |  ____|  ____|   /\    / ____|         (_)     | |    //
+//  | |__  | |__     /  \  | (___   ___ ____ _ ____ | |_   //
+//  |  __| |  __|   / /\ \  \___ \ / __|  __| |  _ \| __|  //
+//  | |    | |____ / ____ \ ____) | (__| |  | | |_) | |    //
+//  |_|    |______/_/    \_\_____/ \___|_|  |_|  __/| |    //
+//                                            | |   | |    //
+//                                            |_|   | |_   //
+//       Website: https://feascript.com/             \__|  //
+
+// Example usage of generalFormPDESolver with a simple mesh and coefficients
+import { generalFormPDESolver } from '../../src/solvers/generalFormPDEScript.js';
+
+// Generate a uniform mesh from 0 to 1 with 20 nodes
+const nNodes = 20;
+const mesh = Array.from({ length: nNodes }, (_, i) => i / (nNodes - 1));
+
+// Coefficient functions for Given equation: d²u/dx² + 10 du/dx = -10 * exp(-200 * (x - 0.5)²), for 0 < x < 1
+const A = x => 1; // Diffusion coefficient
+const B = x => 10; // first derivative term
+const C = x => 0; // No reaction term
+const D = x => -10 * Math.exp(-200 * Math.pow(x - 0.5, 2)); // Source function D(x)
+
+// Dirichlet left boundary conditions: u(0) = 1, and Neumann right boundry conditions: u(1) = 0
+const boundary = {
+    left: { type: 'dirichlet', value: 1 },
+    right: { type: 'neumann', value: 0 }
+};
+
+//Solve
+const u = generalFormPDESolver({ A, B, C, D, mesh, boundary });
+console.log('Mesh nodes:', mesh);
+console.log('Solution u:', u);

examples/generalFormPDEScript/readme.md

@@ -0,0 +1,58 @@
+# Example: 1D General PDE Solver
+
+This example demonstrates how to use the FEAScript 1D general PDE solver for equations of the form:
+
+$$
+\frac{d}{dx}\left(A(x)\frac{du}{dx}\right) + B(x)\frac{du}{dx} + C(x)u = D(x)
+$$
+
+where $u$ is the unknown, $x$ is the position, and $A$, $B$, $C$, $D$ are user-provided coefficient functions.
+
+## How to Run
+
+1. **Set up the mesh:**
+   - Here, we use a uniform mesh from 0 to 1 with 20 nodes.
+2. **Define coefficient functions:**
+   - For the Given equation $-\frac{d^2u}{dx^2} = 1$, set $A(x) = 1$, $B(x) = 10$,
+    $C(x) = 0$, $D(x) = -10 * exp(-200 * (x - 0.5)²)$, for 0 < x < 1.
+3. **Set boundary conditions:**
+   - Dirichlet boundaries: $u(0) = 1$, $u(1) = 0$.
+4. **Solve and print results:**
+   - The solution vector $u$ will be printed for each mesh node.
+
+## Example Code
+
+```js
+const { generalFormPDESolver } = require("../../src/solvers/generalFormPDEScript");
+
+// Generate a uniform mesh from 0 to 1 with 20 nodes
+const nNodes = 20;
+const mesh = Array.from({ length: nNodes }, (_, i) => i / (nNodes - 1));
+
+// Coefficient functions for Given equation: d²u/dx² + 10 du/dx = -10 * exp(-200 * (x - 0.5)²), for 0 < x < 1
+const A = x => 1; // Diffusion coefficient
+const B = x => 0; // No first derivative term
+const C = x => 0; // No reaction term
+const D = x => 1; // Source function D(X)
+
+
+// Dirichlet left boundary conditions: u(0) = 1, and Neumann right boundry conditions: u(1) = 0
+const boundary = {
+    left: { type: 'dirichlet', value: 1 },
+    right: { type: 'neumann', value: 0 }
+};
+
+// Solve
+const u = generalFormPDESolver({ A, B, C, D, mesh, boundary });
+
+console.log('Mesh nodes:', mesh);
+console.log('Solution u:', u);
+
+
+## Output
+
+The script prints the mesh nodes and the solution $u$ at each node. You can modify the coefficient functions and boundary conditions to solve other problems.
+
+---
+
+For more details, see the main project README or documentation.

src/solvers/generalFormPDEScript.js

@@ -0,0 +1,165 @@
+//   ______ ______           _____           _       _     //
+//  |  ____|  ____|   /\    / ____|         (_)     | |    //
+//  | |__  | |__     /  \  | (___   ___ ____ _ ____ | |_   //
+//  |  __| |  __|   / /\ \  \___ \ / __|  __| |  _ \| __|  //
+//  | |    | |____ / ____ \ ____) | (__| |  | | |_) | |    //
+//  |_|    |______/_/    \_\_____/ \___|_|  |_|  __/| |    //
+//                                            | |   | |    //
+//                                            |_|   | |_   //
+//       Website: https://feascript.com/             \__|  //
+
+// generalFormPDEScript.js
+// Solver for general 1D linear differential equations in weak form
+// User provides coefficients A(x), B(x), C(x), D(x) and boundary conditions
+
+/**
+ * General 1D PDE Solver (Weak Form)
+ * Solves equations of the form:
+ * -->  d/dx(A(x) du/dx) + B(x) du/dx + C(x) u = D(x)
+ * using the finite element method (FEM) and user-supplied coefficients.
+ * @param {Object} options - Solver options
+ * @param {function} options.A - Coefficient function A(x) (diffusion)
+ * @param {function} options.B - Coefficient function B(x) (advection)
+ * @param {function} options.C - Coefficient function C(x) (reaction)
+ * @param {function} options.D - Source function D(x)
+ * @param {Array<number>} options.mesh - Mesh nodes (1D array)
+ * @param {Object} options.boundary - Boundary conditions
+ *   { left: {type: 'dirichlet'|'neumann'|'robin', value}, right: {type, value} }
+ *     - Dirichlet: value is the prescribed u
+ *     - Neumann: value is the prescribed flux (A du/dx)
+ *     - Robin: value is { a, b, c } for a*u + b*du/dx = c
+ * @returns {Array<number>} Solution vector u at mesh nodes
+ */
+export function generalFormPDESolver({ A, B, C, D, mesh, boundary }) {
+    // Number of nodes and elements
+    const nNodes = mesh.length;
+    const nElements = nNodes - 1;
+
+    // Initialize global stiffness matrix and load vector
+    const K = Array.from({ length: nNodes }, () => Array(nNodes).fill(0));
+    const F = Array(nNodes).fill(0);
+
+    // Linear basis functions for 1D elements
+    // For each element [x0, x1]
+    for (let e = 0; e < nElements; e++) {
+        const x0 = mesh[e];
+        const x1 = mesh[e + 1];
+        const h = x1 - x0;
+        // Local stiffness matrix and load vector
+        let Ke = [ [0, 0], [0, 0] ];
+        let Fe = [0, 0];
+        // 2-point Gauss quadrature for integration
+        const gaussPoints = [ -1/Math.sqrt(3), 1/Math.sqrt(3) ];
+        const gaussWeights = [ 1, 1 ];
+        for (let gp = 0; gp < 2; gp++) {
+            // Map reference [-1,1] to [x0,x1]
+            const xi = gaussPoints[gp];
+            const w = gaussWeights[gp];
+            const x = ((1 - xi) * x0 + (1 + xi) * x1) / 2;
+            const dx_dxi = h / 2;
+            // Basis functions and derivatives
+            const N = [ (1 - xi) / 2, (1 + xi) / 2 ];
+            const dN_dxi = [ -0.5, 0.5 ];
+            const dN_dx = dN_dxi.map(d => d / dx_dxi);
+            // Coefficients at integration point
+            const a = A(x);
+            const b = B(x);
+            const c = C(x);
+            const d = D(x);
+            // Element matrix assembly (weak form)
+            for (let i = 0; i < 2; i++) {
+                for (let j = 0; j < 2; j++) {
+                    Ke[i][j] += (
+                        a * dN_dx[i] * dN_dx[j] +
+                        b * dN_dx[j] * N[i] -
+                        c * N[i] * N[j]
+                    ) * w * dx_dxi;
+                }
+                Fe[i] += d * N[i] * w * dx_dxi;
+            }
+        }
+        // Assemble into global matrix/vector
+        K[e][e]     += Ke[0][0];
+        K[e][e+1]   += Ke[0][1];
+        K[e+1][e]   += Ke[1][0];
+        K[e+1][e+1] += Ke[1][1];
+        F[e]     += Fe[0];
+        F[e+1]   += Fe[1];
+    }
+
+    // Apply boundary conditions (Dirichlet, Neumann, Robin)
+    // Left boundary
+    if (boundary.left) {
+        if (boundary.left.type === 'dirichlet') {
+            for (let j = 0; j < nNodes; j++) {
+                K[0][j] = 0;
+            }
+            K[0][0] = 1;
+            F[0] = boundary.left.value;
+        } else if (boundary.left.type === 'neumann') {
+            // Add Neumann flux to load vector
+            F[0] += boundary.left.value;
+        } else if (boundary.left.type === 'robin') {
+            // Robin: a*u + b*du/dx = c
+            const { a, b, c } = boundary.left.value;
+            K[0][0] += a;
+            F[0] += c;
+            K[0][0] += b;
+        }
+    }
+    // Right boundary
+    if (boundary.right) {
+        if (boundary.right.type === 'dirichlet') {
+            for (let j = 0; j < nNodes; j++) {
+                K[nNodes-1][j] = 0;
+            }
+            K[nNodes-1][nNodes-1] = 1;
+            F[nNodes-1] = boundary.right.value;
+        } else if (boundary.right.type === 'neumann') {
+            // Add Neumann flux to load vector
+            F[nNodes-1] += boundary.right.value;
+        } else if (boundary.right.type === 'robin') {
+            // Robin: a*u + b*du/dx = c
+            const { a, b, c } = boundary.right.value;
+            K[nNodes-1][nNodes-1] += a;
+            F[nNodes-1] += c;
+            K[nNodes-1][nNodes-1] += b;
+        }
+    }
+
+    // Solve linear system Ku = F (simple Gauss elimination for small systems)
+    function solveLinearSystem(A, b) {
+        // Naive Gauss elimination (for small systems)
+        const n = b.length;
+        const x = Array(n).fill(0);
+        const M = A.map(row => row.slice());
+        const f = b.slice();
+        // Forward elimination
+        for (let i = 0; i < n; i++) {
+            let maxRow = i;
+            for (let k = i+1; k < n; k++) {
+                if (Math.abs(M[k][i]) > Math.abs(M[maxRow][i])) maxRow = k;
+            }
+            [M[i], M[maxRow]] = [M[maxRow], M[i]];
+            [f[i], f[maxRow]] = [f[maxRow], f[i]];
+            for (let k = i+1; k < n; k++) {
+                const c = M[k][i] / M[i][i];
+                for (let j = i; j < n; j++) {
+                    M[k][j] -= c * M[i][j];
+                }
+                f[k] -= c * f[i];
+            }
+        }
+        // Back substitution
+        for (let i = n-1; i >= 0; i--) {
+            let sum = 0;
+            for (let j = i+1; j < n; j++) {
+                sum += M[i][j] * x[j];
+            }
+            x[i] = (f[i] - sum) / M[i][i];
+        }
+        return x;
+    }
+    const u = solveLinearSystem(K, F);
+    return u;
+}