basisFunctions.js

/**
 * ════════════════════════════════════════════════════════════════
 *  FEAScript Core Library
 *  Finite Element Simulation in JavaScript
 *  Version: 0.3.0 | https://feascript.com
 *  MIT License © 2023–2026 FEAScript
 * ════════════════════════════════════════════════════════════════
 */

// Internal imports
import { errorLog } from "../utilities/logging.js";

/**
 * Local node indices on each side of the reference triangle, following the quad side numbering
 * (side 0: bottom, eta = 0; side 1: left, ksi = 0; side 2: hypotenuse, ksi + eta = 1)
 * Nodes are listed in order along each side
 */
export const triangleBoundarySides = {
  linear: {
    0: [0, 2], // Nodes at the bottom side of the reference triangle
    1: [0, 1], // Nodes at the left side of the reference triangle
    2: [1, 2], // Nodes at the hypotenuse of the reference triangle
  },
  quadratic: {
    0: [0, 3, 5], // Nodes at the bottom side of the reference triangle
    1: [0, 1, 2], // Nodes at the left side of the reference triangle
    2: [2, 4, 5], // Nodes at the hypotenuse of the reference triangle
  },
};

/**
 * Local node indices on each side of the 1D reference element
 */
const lineBoundarySides = {
  linear: {
    0: [0], // Node at the left side of the reference element
    1: [1], // Node at the right side of the reference element
  },
  quadratic: {
    0: [0], // Node at the left side of the reference element
    1: [2], // Node at the right side of the reference element
  },
};

/**
 * Local node indices on each side of the reference quadrilateral
 * (side 0: bottom, eta = 0; side 1: left, ksi = 0; side 2: top, eta = 1; side 3: right, ksi = 1)
 * Nodes are listed in order along each side
 */
const quadBoundarySides = {
  linear: {
    0: [0, 2], // Nodes at the bottom side of the reference element
    1: [0, 1], // Nodes at the left side of the reference element
    2: [1, 3], // Nodes at the top side of the reference element
    3: [2, 3], // Nodes at the right side of the reference element
  },
  quadratic: {
    0: [0, 3, 6], // Nodes at the bottom side of the reference element
    1: [0, 1, 2], // Nodes at the left side of the reference element
    2: [2, 5, 8], // Nodes at the top side of the reference element
    3: [6, 7, 8], // Nodes at the right side of the reference element
  },
};

/**
 * Function to get the local node indices on each side of the reference element
 * @param {string} meshDimension - The dimension of the mesh, either '1D' or '2D'
 * @param {string} elementOrder - The order of elements, either 'linear' or 'quadratic'
 * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle'
 * @returns {object} Map from side index to the local node indices on that side
 */
export function getBoundarySides(meshDimension, elementOrder, elementShape = "quad") {
  if (meshDimension === "1D") return lineBoundarySides[elementOrder];
  if (elementShape === "triangle") return triangleBoundarySides[elementOrder];
  return quadBoundarySides[elementOrder];
}

/**
 * Class to handle basis functions and their derivatives based on element configuration
 */
export class BasisFunctions {
  /**
   * Constructor to initialize the BasisFunctions class
   * @param {string} meshDimension - The dimension of the mesh
   * @param {string} [elementOrder='linear'] - The order of elements
   * @param {string} [elementShape='quad'] - The 2D element shape, either 'quad' or 'triangle'
   */
  constructor({ meshDimension, elementOrder = "linear", elementShape = "quad" }) {
    this.meshDimension = meshDimension;
    this.elementOrder = elementOrder;
    this.elementShape = elementShape;
  }

  /**
   * Function to calculate basis functions and their derivatives based on the dimension and order
   * @param {number} ksi - Natural coordinate (for both 1D and 2D), in [0, 1]
   * @param {number} [eta] - Second natural coordinate (only for 2D elements)
   * @param {number} [elementLength] - Physical element length, only required for 1D 'hermiteCubic' elements
   * @returns {object} An object containing:
   *  - basisFunction: Array of evaluated basis functions
   *  - basisFunctionDerivKsi: Array of derivatives of basis functions with respect to ksi
   *  - basisFunctionDerivEta: Array of derivatives of basis functions with respect to eta (only for 2D elements)
   *  - basisFunctionDerivKsi2: Array of second derivatives of basis functions with respect to ksi
   *    (only for 1D 'hermiteCubic' elements)
   *
   * 'hermiteCubic' is a general-purpose interpolation type, not specific to beams: it makes both
   * the value and the slope continuous across elements (unlike the C0 Lagrange types above), which
   * any fourth-order problem needs (e.g. beam or plate bending). eulerBernoulliBeam.js is currently
   * the only place in this codebase that uses it.
   */
  getBasisFunctions(ksi, eta = null, elementLength = null) {
    let basisFunction = [];
    let basisFunctionDerivKsi = [];
    let basisFunctionDerivEta = [];
    let basisFunctionDerivKsi2 = [];

    if (this.meshDimension === "1D") {
      if (this.elementOrder === "linear") {
        // Linear basis functions for 1D elements
        basisFunction[0] = 1 - ksi;
        basisFunction[1] = ksi;

        // Derivatives of basis functions with respect to ksi
        basisFunctionDerivKsi[0] = -1;
        basisFunctionDerivKsi[1] = 1;
      } else if (this.elementOrder === "quadratic") {
        // Quadratic basis functions for 1D elements
        basisFunction[0] = 1 - 3 * ksi + 2 * ksi ** 2;
        basisFunction[1] = 4 * ksi - 4 * ksi ** 2;
        basisFunction[2] = -ksi + 2 * ksi ** 2;

        // Derivatives of basis functions with respect to ksi
        basisFunctionDerivKsi[0] = -3 + 4 * ksi;
        basisFunctionDerivKsi[1] = 4 - 8 * ksi;
        basisFunctionDerivKsi[2] = -1 + 4 * ksi;
      } else if (this.elementOrder === "hermiteCubic") {
        // Cubic Hermite basis functions for 1D Euler-Bernoulli beam elements
        // DOF order: [w1, theta1, w2, theta2], with w = transverse deflection and
        // theta = dw/dx the (physical) slope. The h-scaling on the theta-associated
        // functions converts the interpolated dw/dksi at the nodes into the actual
        // rotation DOF, so no extra scaling is required when mapping derivatives to x
        if (elementLength === null) {
          errorLog("elementLength is required to evaluate 'hermiteCubic' basis functions");
          return;
        }
        const h = elementLength;

        basisFunction[0] = 2 * ksi ** 3 - 3 * ksi ** 2 + 1;
        basisFunction[1] = h * (ksi ** 3 - 2 * ksi ** 2 + ksi);
        basisFunction[2] = -2 * ksi ** 3 + 3 * ksi ** 2;
        basisFunction[3] = h * (ksi ** 3 - ksi ** 2);

        // First derivatives of basis functions with respect to ksi
        basisFunctionDerivKsi[0] = 6 * ksi ** 2 - 6 * ksi;
        basisFunctionDerivKsi[1] = h * (3 * ksi ** 2 - 4 * ksi + 1);
        basisFunctionDerivKsi[2] = -6 * ksi ** 2 + 6 * ksi;
        basisFunctionDerivKsi[3] = h * (3 * ksi ** 2 - 2 * ksi);

        // Second derivatives of basis functions with respect to ksi
        basisFunctionDerivKsi2[0] = 12 * ksi - 6;
        basisFunctionDerivKsi2[1] = h * (6 * ksi - 4);
        basisFunctionDerivKsi2[2] = -12 * ksi + 6;
        basisFunctionDerivKsi2[3] = h * (6 * ksi - 2);
      }
    } else if (this.meshDimension === "2D") {
      if (eta === null) {
        errorLog("Eta coordinate is required for 2D elements");
        return;
      }

      if (this.elementShape === "triangle") {
        // Reference triangle with vertices (0,0), (0,1), (1,0); area coordinate l1 = 1 - ksi - eta
        const l1 = 1 - ksi - eta;

        if (this.elementOrder === "linear") {
          // Nodes numbered from bottom to top, then left to right (same rule as the quads):
          // 1
          // | \
          // |   \
          // 0 --- 2
          basisFunction[0] = l1;
          basisFunction[1] = eta;
          basisFunction[2] = ksi;

          basisFunctionDerivKsi[0] = -1;
          basisFunctionDerivKsi[1] = 0;
          basisFunctionDerivKsi[2] = 1;

          basisFunctionDerivEta[0] = -1;
          basisFunctionDerivEta[1] = 1;
          basisFunctionDerivEta[2] = 0;
        } else if (this.elementOrder === "quadratic") {
          // Nodes numbered from bottom to top, then left to right (same rule as the quads):
          // 2
          // | \
          // 1   4
          // |     \
          // 0--3--5
          basisFunction[0] = l1 * (2 * l1 - 1);
          basisFunction[1] = 4 * l1 * eta;
          basisFunction[2] = eta * (2 * eta - 1);
          basisFunction[3] = 4 * l1 * ksi;
          basisFunction[4] = 4 * ksi * eta;
          basisFunction[5] = ksi * (2 * ksi - 1);

          basisFunctionDerivKsi[0] = 1 - 4 * l1;
          basisFunctionDerivKsi[1] = -4 * eta;
          basisFunctionDerivKsi[2] = 0;
          basisFunctionDerivKsi[3] = 4 * (l1 - ksi);
          basisFunctionDerivKsi[4] = 4 * eta;
          basisFunctionDerivKsi[5] = 4 * ksi - 1;

          basisFunctionDerivEta[0] = 1 - 4 * l1;
          basisFunctionDerivEta[1] = 4 * (l1 - eta);
          basisFunctionDerivEta[2] = 4 * eta - 1;
          basisFunctionDerivEta[3] = -4 * ksi;
          basisFunctionDerivEta[4] = 4 * ksi;
          basisFunctionDerivEta[5] = 0;
        }
      } else if (this.elementOrder === "linear") {
        // Linear basis functions for 2D elements
        function l1(c) {
          return 1 - c;
        }
        function l2(c) {
          return c;
        }
        function dl1() {
          return -1;
        }
        function dl2() {
          return 1;
        }

        // Evaluate basis functions at (ksi, eta)
        basisFunction[0] = l1(ksi) * l1(eta);
        basisFunction[1] = l1(ksi) * l2(eta);
        basisFunction[2] = l2(ksi) * l1(eta);
        basisFunction[3] = l2(ksi) * l2(eta);

        // Derivatives with respect to ksi
        basisFunctionDerivKsi[0] = dl1() * l1(eta);
        basisFunctionDerivKsi[1] = dl1() * l2(eta);
        basisFunctionDerivKsi[2] = dl2() * l1(eta);
        basisFunctionDerivKsi[3] = dl2() * l2(eta);

        // Derivatives with respect to eta
        basisFunctionDerivEta[0] = l1(ksi) * dl1();
        basisFunctionDerivEta[1] = l1(ksi) * dl2();
        basisFunctionDerivEta[2] = l2(ksi) * dl1();
        basisFunctionDerivEta[3] = l2(ksi) * dl2();
      } else if (this.elementOrder === "quadratic") {
        // Quadratic basis functions for 2D elements
        function l1(c) {
          return 2 * c ** 2 - 3 * c + 1;
        }
        function l2(c) {
          return -4 * c ** 2 + 4 * c;
        }
        function l3(c) {
          return 2 * c ** 2 - c;
        }
        function dl1(c) {
          return 4 * c - 3;
        }
        function dl2(c) {
          return -8 * c + 4;
        }
        function dl3(c) {
          return 4 * c - 1;
        }

        // Evaluate basis functions at (ksi, eta)
        basisFunction[0] = l1(ksi) * l1(eta);
        basisFunction[1] = l1(ksi) * l2(eta);
        basisFunction[2] = l1(ksi) * l3(eta);
        basisFunction[3] = l2(ksi) * l1(eta);
        basisFunction[4] = l2(ksi) * l2(eta);
        basisFunction[5] = l2(ksi) * l3(eta);
        basisFunction[6] = l3(ksi) * l1(eta);
        basisFunction[7] = l3(ksi) * l2(eta);
        basisFunction[8] = l3(ksi) * l3(eta);

        // Derivatives with respect to ksi
        basisFunctionDerivKsi[0] = dl1(ksi) * l1(eta);
        basisFunctionDerivKsi[1] = dl1(ksi) * l2(eta);
        basisFunctionDerivKsi[2] = dl1(ksi) * l3(eta);
        basisFunctionDerivKsi[3] = dl2(ksi) * l1(eta);
        basisFunctionDerivKsi[4] = dl2(ksi) * l2(eta);
        basisFunctionDerivKsi[5] = dl2(ksi) * l3(eta);
        basisFunctionDerivKsi[6] = dl3(ksi) * l1(eta);
        basisFunctionDerivKsi[7] = dl3(ksi) * l2(eta);
        basisFunctionDerivKsi[8] = dl3(ksi) * l3(eta);

        // Derivatives with respect to eta
        basisFunctionDerivEta[0] = l1(ksi) * dl1(eta);
        basisFunctionDerivEta[1] = l1(ksi) * dl2(eta);
        basisFunctionDerivEta[2] = l1(ksi) * dl3(eta);
        basisFunctionDerivEta[3] = l2(ksi) * dl1(eta);
        basisFunctionDerivEta[4] = l2(ksi) * dl2(eta);
        basisFunctionDerivEta[5] = l2(ksi) * dl3(eta);
        basisFunctionDerivEta[6] = l3(ksi) * dl1(eta);
        basisFunctionDerivEta[7] = l3(ksi) * dl2(eta);
        basisFunctionDerivEta[8] = l3(ksi) * dl3(eta);
      }
    }

    return { basisFunction, basisFunctionDerivKsi, basisFunctionDerivEta, basisFunctionDerivKsi2 };
  }
}