Finite Element Method — Structural Analysis from First Principles

1. Linear Elasticity and Governing Equations

Linear elasticity assumes small strains, linear material response (Hooke's law), and no plastic deformation. The equilibrium PDE in 3D is:

2. Weak Form and Galerkin Method

3. Shape Functions and Isoparametric Elements

4. Element Stiffness Matrix K_e

5. Global Assembly and Boundary Conditions

6. Stress Recovery and Von Mises

7. JavaScript CST Element Solver

// Minimal 2D plane-stress FEM with CST triangles
class FEM2D {
  constructor(nodes, elements, E, nu, thickness = 1) {
    this.nodes    = nodes;    // [[x,y], ...]
    this.elements = elements; // [[n0,n1,n2], ...]
    this.t = thickness;
    this.D = planeStressD(E, nu);
    const ndof = 2 * nodes.length;
    this.K = Array.from({length: ndof}, () => new Float64Array(ndof));
    this.f = new Float64Array(ndof);
  }

  assemble() {
    for (const [n0, n1, n2] of this.elements) {
      const [Ke, Be, A] = cstElement(this.nodes[n0], this.nodes[n1],
                                         this.nodes[n2], this.D, this.t);
      const dofs = [2*n0, 2*n0+1, 2*n1, 2*n1+1, 2*n2, 2*n2+1];
      for (let i = 0; i < 6; i++)
        for (let j = 0; j < 6; j++)
          this.K[dofs[i]][dofs[j]] += Ke[i][j];
    }
  }

  applyFixed(nodeIndex, axis) { // axis: 0=x, 1=y
    const d = 2 * nodeIndex + axis;
    const n = this.K.length;
    for (let j = 0; j < n; j++) this.K[d][j] = this.K[j][d] = 0;
    this.K[d][d] = 1;
    this.f[d] = 0;
  }

  applyForce(nodeIndex, fx, fy) {
    this.f[2*nodeIndex]   += fx;
    this.f[2*nodeIndex+1] += fy;
  }

  solve() { // Naive Gaussian elimination (use sparse solver in production)
    const A = this.K.map(r => Array.from(r));
    const b = Array.from(this.f);
    const n = b.length;
    for (let p = 0; p < n; p++) {
      for (let i = p+1; i < n; i++) {
        const f = A[i][p] / A[p][p];
        for (let j = p; j < n; j++) A[i][j] -= f*A[p][j];
        b[i] -= f*b[p];
      }
    }
    const u = new Float64Array(n);
    for (let i = n-1; i >= 0; i--) {
      u[i] = b[i];
      for (let j = i+1; j < n; j++) u[i] -= A[i][j]*u[j];
      u[i] /= A[i][i];
    }
    return u;
  }
}

function planeStressD(E, nu) {
  const c = E / (1 - nu*nu);
  return [[c, c*nu, 0], [c*nu, c, 0], [0, 0, c*(1-nu)/2]];
}

function cstElement([x1,y1], [x2,y2], [x3,y3], D, t) {
  const A = 0.5 * Math.abs((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1));
  const b = [y2-y3, y3-y1, y1-y2];
  const c = [x3-x2, x1-x3, x2-x1];
  // B is 3×6
  const B = [
    [b[0],0,  b[1],0,  b[2],0  ],
    [0,  c[0],0,  c[1],0,  c[2]],
    [c[0],b[0],c[1],b[1],c[2],b[2]]
  ].map(row => row.map(v => v / (2*A)));
  // Ke = B^T D B A t  (6×6)
  const DB = matMul(D, B);  // 3×6
  const Ke = matMul(transpose(B), DB).map(row => row.map(v => v*A*t));
  return [Ke, B, A];
}

8. Element Types and Software

Commercial FEM codes — Ansys, Abaqus, NASTRAN, LS-DYNA — add contact mechanics, nonlinear materials, geometric nonlinearity, dynamic (modal and transient) analysis, and multiphysics coupling. Open-source alternatives include FEniCS (Python, FEniCS forms language) and Code_Aster (general) and deal.II (C++ templates for high-order adaptive methods).