// Triangular mesh of the projected disk, as in cymsim/mesh.py: nodes on concentric rings,
// annuli stitched by an angular merge, optional Laplacian smoothing, and the constant
// gradient operators of the linear (CST) interpolation on every triangle.

export interface Mesh {
  n: number; e: number;
  xy: Float64Array;        // (N, 2)
  tri: Uint32Array;        // (E, 3) counter-clockwise
  boundary: Uint32Array;   // rim node indices in angular order
  area: Float64Array;      // (E,)
  Gx: Float64Array;        // (E, 3): d phi/dx = sum_i Gx[e, i] phi_i
  Gy: Float64Array;
  centroid: Float64Array;  // (E, 2)
}

function stitch(aIdx: Uint32Array, aAng: Float64Array, bIdx: Uint32Array, bAng: Float64Array, out: number[]) {
  const na = aIdx.length, nb = bIdx.length, twoPi = 2 * Math.PI;
  const bAngU = Array.from(bAng, (t) => (((t - aAng[0]) % twoPi) + twoPi) % twoPi + aAng[0]);
  const order = Array.from(bAngU.keys()).sort((i, j) => bAngU[i] - bAngU[j]);
  const bi = order.map((i) => bIdx[i]), ba = order.map((i) => bAngU[i]);
  const ta = Array.from(aAng); ta.push(aAng[0] + twoPi);
  const tb = ba.slice(); tb.push(ba[0] + twoPi);
  let i = 0, j = 0;
  while (i < na || j < nb) {
    const nextA = i < na ? ta[i + 1] : Infinity, nextB = j < nb ? tb[j + 1] : Infinity;
    if (nextA <= nextB) { out.push(aIdx[i % na], bi[j % nb], aIdx[(i + 1) % na]); i++; }
    else { out.push(aIdx[i % na], bi[j % nb], bi[(j + 1) % nb]); j++; }
  }
}

function smooth(xy: Float64Array, tri: Uint32Array, fixed: Uint32Array, nIter: number): Float64Array {
  const N = xy.length / 2, E = tri.length / 3;
  const src: number[] = [], dst: number[] = [];
  for (let e = 0; e < E; e++) {
    for (let k = 0; k < 3; k++) {
      const a = tri[3 * e + k], b = tri[3 * e + ((k + 1) % 3)];
      src.push(a, b); dst.push(b, a);
    }
  }
  const deg = new Float64Array(N);
  for (const s of src) deg[s]++;
  const free = new Uint8Array(N).fill(1);
  for (const f of fixed) free[f] = 0;
  let cur = Float64Array.from(xy);
  for (let it = 0; it < nIter; it++) {
    const acc = new Float64Array(2 * N);
    for (let k = 0; k < src.length; k++) { acc[2 * src[k]] += cur[2 * dst[k]]; acc[2 * src[k] + 1] += cur[2 * dst[k] + 1]; }
    const next = Float64Array.from(cur);
    for (let n = 0; n < N; n++) if (free[n]) { next[2 * n] = acc[2 * n] / deg[n]; next[2 * n + 1] = acc[2 * n + 1] / deg[n]; }
    cur = next;
  }
  return cur;
}

export function diskMesh(radius: number, elementSize: number, smoothing = 1): Mesh {
  const nRings = Math.max(2, Math.round(radius / elementSize));
  const dr = radius / nRings;
  const nodes: number[] = [0, 0];
  let prevIdx = new Uint32Array([0]), prevAng = new Float64Array([0]);
  const tris: number[] = [];
  for (let k = 1; k <= nRings; k++) {
    const r = k * dr;
    const nk = Math.max(6, Math.round((2 * Math.PI * r) / dr));
    const ang = new Float64Array(nk), idx = new Uint32Array(nk);
    const base = nodes.length / 2;
    for (let j = 0; j < nk; j++) {
      ang[j] = (2 * Math.PI * (j + 0.5 * (k % 2))) / nk;
      idx[j] = base + j;
      nodes.push(r * Math.cos(ang[j]), r * Math.sin(ang[j]));
    }
    if (k === 1) for (let j = 0; j < nk; j++) tris.push(0, idx[j], idx[(j + 1) % nk]);
    else stitch(prevIdx, prevAng, idx, ang, tris);
    prevIdx = idx; prevAng = ang;
  }
  let xy: Float64Array = Float64Array.from(nodes);
  const tri = Uint32Array.from(tris);
  if (smoothing > 0) xy = smooth(xy, tri, prevIdx, smoothing);
  const E = tri.length / 3;
  for (let e = 0; e < E; e++) {
    const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
    const s = (xy[2 * b] - xy[2 * a]) * (xy[2 * c + 1] - xy[2 * a + 1]) - (xy[2 * c] - xy[2 * a]) * (xy[2 * b + 1] - xy[2 * a + 1]);
    if (s < 0) { tri[3 * e + 1] = c; tri[3 * e + 2] = b; }
  }
  return makeMesh(xy, tri, prevIdx);
}

export function makeMesh(xy: Float64Array, tri: Uint32Array, boundary: Uint32Array): Mesh {
  const N = xy.length / 2, E = tri.length / 3;
  const area = new Float64Array(E), Gx = new Float64Array(3 * E), Gy = new Float64Array(3 * E), centroid = new Float64Array(2 * E);
  for (let e = 0; e < E; e++) {
    const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
    const x0 = xy[2 * a], y0 = xy[2 * a + 1], x1 = xy[2 * b], y1 = xy[2 * b + 1], x2 = xy[2 * c], y2 = xy[2 * c + 1];
    const A = 0.5 * ((x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0));
    if (!(A > 0)) throw new Error(`element ${e} has non-positive area`);
    area[e] = A;
    Gx[3 * e] = (y1 - y2) / (2 * A); Gx[3 * e + 1] = (y2 - y0) / (2 * A); Gx[3 * e + 2] = (y0 - y1) / (2 * A);
    Gy[3 * e] = (x2 - x1) / (2 * A); Gy[3 * e + 1] = (x0 - x2) / (2 * A); Gy[3 * e + 2] = (x1 - x0) / (2 * A);
    centroid[2 * e] = (x0 + x1 + x2) / 3; centroid[2 * e + 1] = (y0 + y1 + y2) / 3;
  }
  return { n: N, e: E, xy, tri, boundary, area, Gx, Gy, centroid };
}

/** Containing element and barycentric weights of a point; throws when outside the mesh. */
export function locate(m: Mesh, px: number, py: number): { element: number; nodes: [number, number, number]; weights: [number, number, number] } {
  const { xy, tri, area } = m;
  for (let e = 0; e < m.e; e++) {
    const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
    const x0 = xy[2 * a], y0 = xy[2 * a + 1], x1 = xy[2 * b], y1 = xy[2 * b + 1], x2 = xy[2 * c], y2 = xy[2 * c + 1];
    const a2 = 2 * area[e];
    const l1 = ((y1 - y2) * (px - x2) + (x2 - x1) * (py - y2)) / a2;
    const l2 = ((y2 - y0) * (px - x2) + (x0 - x2) * (py - y2)) / a2;
    const l3 = 1 - l1 - l2;
    if (l1 >= -1e-9 && l2 >= -1e-9 && l3 >= -1e-9) return { element: e, nodes: [a, b, c], weights: [l1, l2, l3] };
  }
  throw new Error(`point (${px.toFixed(4)}, ${py.toFixed(4)}) lies outside the mesh`);
}

export function edgeLengths(m: Mesh): Float64Array {
  const L = new Float64Array(3 * m.e);
  for (let e = 0; e < m.e; e++) {
    for (let k = 0; k < 3; k++) {
      const a = m.tri[3 * e + k], b = m.tri[3 * e + ((k + 1) % 3)];
      L[3 * e + k] = Math.hypot(m.xy[2 * a] - m.xy[2 * b], m.xy[2 * a + 1] - m.xy[2 * b + 1]);
    }
  }
  return L;
}