Port of https://github.com/matthias-research/pages/blob/master/tenMinutePhysics/BlenderTetPlugin.py to TypeScript

blender-tetrahedron-generation-plugin.ts

import { MeshBVH, acceleratedRaycast } from 'three-mesh-bvh'
import * as THREE from 'three';

THREE.Mesh.prototype.raycast = acceleratedRaycast;

const raycaster = new THREE.Raycaster();
function print(...args: any[]) {
  console.log(...args);
}
function len<T>(a: ArrayLike<T>) {
  return a.length;
}
function int(val: number) {
  return 0 | val;
}
function random() {
  return Math.random();
}
class mathutils {
  static Vector(v: number[]) {
    return new THREE.Vector3(...v);
  }
}
function max(a: number, b: number) {
  return Math.max(a, b);
}

// ------------------------------------------------
const dirs = [
  new THREE.Vector3(1.0, 0.0, 0.0),
  new THREE.Vector3(0.0, -1.0, 0.0),
  new THREE.Vector3(0.0, 1.0, 0.0),
  new THREE.Vector3(0.0, -1.0, 0.0),
  new THREE.Vector3(0.0, 0.0, 1.0),
  new THREE.Vector3(0.0, 0.0, -1.0)]

type U32x3 = [number, number, number];
const tetFaces: U32x3[] = [[2, 1, 0], [0, 1, 3], [1, 2, 3], [2, 0, 3]]

//

function ray_cast(bvh: MeshBVH, origin: THREE.Vector3, dir: THREE.Vector3): [THREE.Vector3?, THREE.Vector3?, number?, number?] {
  raycaster.ray.origin.copy(origin);
  raycaster.ray.direction.copy(dir);
  const hit: THREE.Intersection = (bvh.raycast(raycaster.ray, THREE.DoubleSide) as Array<THREE.Intersection>).sort((a, b) => a.distance - b.distance)[0] || {};
  return [hit.point, hit.face?.normal, hit.index, hit.distance];
}

const TMP = [
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3(),
  new THREE.Vector3()
]

// ------------------------------------------------

export function isInside(tree: MeshBVH, p: THREE.Vector3, minDist = 0.0) {
  let numIn = 0;
  for (let i = 0; i < 6; i++) {
    const [_location, normal, _index, distance ] = ray_cast(tree, p, dirs[i])
    if (normal) {
      if (normal.dot(dirs[i]) > 0.0) {
        numIn = numIn + 1;
      }
      if (minDist > 0.0 && distance != undefined && distance < minDist) {
        return false;
      }
    }
  }
  return numIn > 3;
};
export function getCircumCenter(p0: THREE.Vector3, p1: THREE.Vector3, p2: THREE.Vector3, p3: THREE.Vector3) {
  const b = TMP[0].subVectors(p1, p0);
  const c = TMP[1].subVectors(p2, p0);
  const d = TMP[2].subVectors(p3, p0);
  const det = 2.0 * ((b.x * (c.y * d.z - c.z * d.y) - b.y * (c.x * d.z - c.z * d.x)) + b.z * (c.x * d.y - c.y * d.x));
  if (det == 0.0) {
    return p0.clone();
  }
  else {
    function crossDot(out: THREE.Vector3, v0: THREE.Vector3, v1: THREE.Vector3, v2: THREE.Vector3) {
      return out.crossVectors(v0, v1).multiplyScalar(v2.dot(v2))
    }
    const v = crossDot(TMP[3], c, d, b).add(crossDot(TMP[4], d, b, c)).add(crossDot(TMP[5], b, c, d))
    v.divideScalar(det);
    return p0.clone().add(v);
  }
};
export function tetQuality(p0: THREE.Vector3, p1: THREE.Vector3, p2: THREE.Vector3, p3: THREE.Vector3) {
  const d0 = TMP[0].subVectors(p1, p0);
  const d1 = TMP[1].subVectors(p2, p0);
  const d2 = TMP[2].subVectors(p3, p0);
  const d3 = TMP[3].subVectors(p2, p1);
  const d4 = TMP[4].subVectors(p3, p2);
  const d5 = TMP[5].subVectors(p1, p3);
  const s0 = d0.length();
  const s1 = d1.length();
  const s2 = d2.length();
  const s3 = d3.length();
  const s4 = d4.length();
  const s5 = d5.length();
  const ms = (s0 * s0 + s1 * s1 + s2 * s2 + s3 * s3 + s4 * s4 + s5 * s5) / 6.0;
  const rms = Math.sqrt(ms);
  const s = 12.0 / Math.sqrt(2.0);
  const vol = d0.dot(d1.cross(d2)) / 6.0;
  return (s * vol) / (rms * rms * rms);
};
export function compareEdges(e0: number[], e1: number[]) {
  if (e0[0] < e1[0] || e0[0] == e1[0] && e0[1] < e1[1]) {
    return -(1);
  }
  else {
    return 1;
  }
};
export function equalEdges(e0: number[], e1: number[]) {
  return e0[0] == e1[0] && e0[1] == e1[1];
};
export function createTetIds(verts: THREE.Vector3[], tree: MeshBVH, minQuality: number) {
  const tetIds: number[] = [];
  const neighbors: number[] = [];
  const tetMarks: number[] = [];
  let tetMark = 0;
  let firstFreeTet = -(1);
  const planesN: THREE.Vector3[] = [];
  const planesD: number[] = [];
  const firstBig = len(verts) - 4;
  tetIds.push(firstBig);
  tetIds.push(firstBig + 1);
  tetIds.push(firstBig + 2);
  tetIds.push(firstBig + 3);
  tetMarks.push(0);
  for (let i = 0; i < 4; i++) {
    neighbors.push(-(1));
    const p0 = verts[firstBig + tetFaces[i][0]];
    const p1 = verts[firstBig + tetFaces[i][1]];
    const p2 = verts[firstBig + tetFaces[i][2]];
    const n = TMP[0].subVectors(p1, p0).cross(TMP[1].subVectors(p2, p0));
    n.normalize();
    planesN.push(n.clone());
    planesD.push(p0.dot(n));
  }
  const center = mathutils.Vector([0.0, 0.0, 0.0]);
  print(' ------------- tetrahedralization ------------------- ');
  for (let i = 0; i < firstBig; i++) {
    const p = verts[i];
    if ((i % 100) == 0) {
      print('inserting vert', i + 1, 'of', firstBig);
    }
    let tetNr = 0;
    while (tetIds[4 * tetNr] < 0) {
      tetNr = tetNr + 1;
    }
    tetMark = tetMark + 1;
    let found = false;
    while (!(found)) {
      if (tetNr < 0 || tetMarks[tetNr] == tetMark) {
        break;
      }
      tetMarks[tetNr] = tetMark;
      const id0 = tetIds[4 * tetNr];
      const id1 = tetIds[4 * tetNr + 1];
      const id2 = tetIds[4 * tetNr + 2];
      const id3 = tetIds[4 * tetNr + 3];
      center.addVectors(verts[id0], verts[id1]).add(verts[id2]).add(verts[id3]).multiplyScalar(0.25);
      let minT = Number.POSITIVE_INFINITY;
      let minFaceNr = -(1);
      for (let j = 0; j < 4; j++) {
        const n = planesN[4 * tetNr + j];
        const d = planesD[4 * tetNr + j];
        const hp = n.dot(p) - d;
        const hc = n.dot(center) - d;
        let t = hp - hc;
        if (t == 0) {
          continue;
        }
        t = -(hc) / t;
        if (t >= 0.0 && t < minT) {
          minT = t;
          minFaceNr = j;
        }
      }
      if (minT >= 1.0) {
        found = true;
      }
      else {
        tetNr = neighbors[4 * tetNr + minFaceNr];
      }
    }
    if (!(found)) {
      print('*********** failed to insert vertex');
      continue;
    }
    tetMark = tetMark + 1;
    const violatingTets: number[] = [];
    const stack = [tetNr];
    while (len(stack) != 0) {
      const tetNr = stack.pop()!;
      if (tetMarks[tetNr] == tetMark) {
        continue;
      }
      tetMarks[tetNr] = tetMark;
      violatingTets.push(tetNr);
      for (let j = 0; j < 4; j++) {
        const n = neighbors[4 * tetNr + j];
        if (n < 0 || tetMarks[n] == tetMark) {
          continue;
        }
        const id0 = tetIds[4 * n];
        const id1 = tetIds[4 * n + 1];
        const id2 = tetIds[4 * n + 2];
        const id3 = tetIds[4 * n + 3];
        const c = getCircumCenter(verts[id0], verts[id1], verts[id2], verts[id3]);
        const r = TMP[0].subVectors(verts[id0], c).length();
        if (TMP[1].subVectors(p, c).length() < r) {
          stack.push(n);
        }
      }
    }
    const edges = [];
    for (let j = 0; j < len(violatingTets); j++) {
      const tetNr = violatingTets[j];
      const ids = [0, 0, 0, 0];
      const ns = [0, 0, 0, 0];
      for (let k = 0; k < 4; k++) {
        ids[k] = tetIds[4 * tetNr + k];
        ns[k] = neighbors[4 * tetNr + k];
      }
      tetIds[4 * tetNr] = -(1);
      tetIds[4 * tetNr + 1] = firstFreeTet;
      firstFreeTet = tetNr;
      for (let k = 0; k < 4; k++) {
        const n = ns[k];
        if (n >= 0 && tetMarks[n] == tetMark) {
          continue;
        }
        let newTetNr = firstFreeTet;
        if (newTetNr >= 0) {
          firstFreeTet = tetIds[4 * firstFreeTet + 1];
        }
        else {
          newTetNr = int(len(tetIds) / 4);
          tetMarks.push(0);
          for (let l = 0; l < 4; l++) {
            tetIds.push(-(1));
            neighbors.push(-(1));
            planesN.push(mathutils.Vector([0.0, 0.0, 0.0]));
            planesD.push(0.0);
          }
        }
        const id0 = ids[tetFaces[k][2]];
        const id1 = ids[tetFaces[k][1]];
        const id2 = ids[tetFaces[k][0]];
        tetIds[4 * newTetNr] = id0;
        tetIds[4 * newTetNr + 1] = id1;
        tetIds[4 * newTetNr + 2] = id2;
        tetIds[4 * newTetNr + 3] = i;
        neighbors[4 * newTetNr] = n;
        if (n >= 0) {
          for (let l = 0; l < 4; l++) {
            if (neighbors[4 * n + l] == tetNr) {
              neighbors[4 * n + l] = newTetNr;
            }
          }
        }
        neighbors[4 * newTetNr + 1] = -(1);
        neighbors[4 * newTetNr + 2] = -(1);
        neighbors[4 * newTetNr + 3] = -(1);
        for (let l = 0; l < 4; l++) {
          const p0 = verts[tetIds[4 * newTetNr + tetFaces[l][0]]];
          const p1 = verts[tetIds[4 * newTetNr + tetFaces[l][1]]];
          const p2 = verts[tetIds[4 * newTetNr + tetFaces[l][2]]];
          const newN = TMP[0].subVectors(p1, p0).cross(TMP[1].subVectors(p2, p0));
          newN.normalize();
          planesN[4 * newTetNr + l] = newN.clone();
          planesD[4 * newTetNr + l] = newN.dot(p0);
        }
        if (id0 < id1) {
          edges.push([id0, id1, newTetNr, 1]);
        }
        else {
          edges.push([id1, id0, newTetNr, 1]);
        }
        if (id1 < id2) {
          edges.push([id1, id2, newTetNr, 2]);
        }
        else {
          edges.push([id2, id1, newTetNr, 2]);
        }
        if (id2 < id0) {
          edges.push([id2, id0, newTetNr, 3]);
        }
        else {
          edges.push([id0, id2, newTetNr, 3]);
        }
      }
    }
    const sortedEdges = edges.sort(compareEdges);
    let nr = 0;
    const numEdges = len(sortedEdges);
    while (nr < numEdges) {
      const e0 = sortedEdges[nr];
      nr = nr + 1;
      if (nr < numEdges && equalEdges(sortedEdges[nr], e0)) {
        const e1 = sortedEdges[nr];
        const id0 = tetIds[4 * e0[2]];
        const id1 = tetIds[4 * e0[2] + 1];
        const id2 = tetIds[4 * e0[2] + 2];
        const id3 = tetIds[4 * e0[2] + 3];
        const jd0 = tetIds[4 * e1[2]];
        const jd1 = tetIds[4 * e1[2] + 1];
        const jd2 = tetIds[4 * e1[2] + 2];
        const jd3 = tetIds[4 * e1[2] + 3];
        neighbors[4 * e0[2] + e0[3]] = e1[2];
        neighbors[4 * e1[2] + e1[3]] = e0[2];
        nr = nr + 1;
      }
    }
  }
  const numTets = int(len(tetIds) / 4);
  let num = 0;
  let numBad = 0;
  for (let i = 0; i < numTets; i++) {
    const id0 = tetIds[4 * i];
    const id1 = tetIds[4 * i + 1];
    const id2 = tetIds[4 * i + 2];
    const id3 = tetIds[4 * i + 3];
    if (id0 < 0 || id0 >= firstBig || id1 >= firstBig || id2 >= firstBig || id3 >= firstBig) {
      continue;
    }
    const p0 = verts[id0];
    const p1 = verts[id1];
    const p2 = verts[id2];
    const p3 = verts[id3];
    const quality = tetQuality(p0, p1, p2, p3);
    if (quality < minQuality) {
      numBad = numBad + 1;
      continue;
    }
    center.addVectors(p0, p1).add(p2).add(p3).multiplyScalar(0.25);
    if (!(isInside(tree, center))) {
      continue;
    }
    tetIds[num] = id0;
    num = num + 1;
    tetIds[num] = id1;
    num = num + 1;
    tetIds[num] = id2;
    num = num + 1;
    tetIds[num] = id3;
    num = num + 1;
  }
  tetIds.splice(num, tetIds.length - num);
  print(numBad, 'bad tets deleted');
  print(int(len(tetIds) / 4), 'tets created');
  return tetIds;
};

export const randEps = function () {
  const eps = 0.0001;
  return -(eps) + (2.0 * random()) * eps;
};

export function createTets(geometry: THREE.BufferGeometry, resolution: number, minQuality: number, oneFacePerTet: boolean, scale = 1.0) {

  const tree = new MeshBVH(geometry);
  const bm = {
    verts: new Array<THREE.Vector3>(),
    faces: new Array<number[]>
  }
  const tetVerts: THREE.Vector3[] = [];
  const vertex = geometry.getAttribute('position').array;
  for (let i = 0; i < vertex.length; i += 3) {
    const v = { co: vertex.subarray(i, i + 3) };
    tetVerts.push(mathutils.Vector([v.co[0] + randEps(), v.co[1] + randEps(), v.co[2] + randEps()]));
  }
  const inf = Number.POSITIVE_INFINITY;
  const center = mathutils.Vector([0.0, 0.0, 0.0]);
  const bmin = [inf, inf, inf];
  const bmax = [-inf, -inf, -inf];
  for (const p of tetVerts) {
    center.add(p);
    const _p = p.toArray();
    for (let i = 0; i < 3; i++) {
      bmin[i] = Math.min(bmin[i], _p[i]);
      bmax[i] = Math.max(bmax[i], _p[i]);
    }
  }
  center.divideScalar(len(tetVerts));
  let radius = 0.0;
  for (const p of tetVerts) {
    const d = TMP[0].subVectors(p, center).length();
    radius = max(radius, d);
  }
  if (resolution > 0) {
    const dims = bmax.map((a, i) => a - bmin[i]);
    const dim = max(dims[0], max(dims[1], dims[2]));
    const h = dim / resolution;
    for (let xi = 0; xi < int(dims[0] / h) + 1; xi++) {
      const x = (bmin[0] + xi * h) + randEps();
      for (let yi = 0; yi < int(dims[1] / h) + 1; yi++) {
        const y = (bmin[1] + yi * h) + randEps();
        for (let zi = 0; zi < int(dims[2] / h) + 1; zi++) {
          const z = (bmin[2] + zi * h) + randEps();
          const p = mathutils.Vector([x, y, z]);
          if (isInside(tree, p, 0.5 * h)) {
            tetVerts.push(p);
          }
        }
      }
    }
  }
  const s = 5.0 * radius;
  tetVerts.push(mathutils.Vector([-s, 0.0, -s]));
  tetVerts.push(mathutils.Vector([s, 0.0, -s]));
  tetVerts.push(mathutils.Vector([0.0, s, s]));
  tetVerts.push(mathutils.Vector([0.0, -s, s]));
  const faces = createTetIds(tetVerts, tree, minQuality);
  const numTets = int(len(faces) / 4);
  if (oneFacePerTet) {
    const numSrcPoints = len(vertex);
    const numPoints = len(tetVerts) - 4;
    for (let i = 0; i < numSrcPoints; i++) {
      const co = vertex.subarray(i * 3, i * 3 + 3);
      bm.verts.push(new THREE.Vector3(...co));
    }
    for (let i = numSrcPoints; i < numPoints; i++) {
      const p = tetVerts[i];
      bm.verts.push(new THREE.Vector3(p.x, p.y, p.z));
    }
  }
  else {
    for (let i = 0; i < numTets; i++) {
      center.addVectors(tetVerts[faces[4 * i]], tetVerts[faces[4 * i + 1]]).add(tetVerts[faces[4 * i + 2]]).add(tetVerts[faces[4 * i + 3]]).multiplyScalar(0.25);
      for (let j = 0; j < 4; j++) {
        for (let k = 0; k < 3; k++) {
          let p = tetVerts[faces[4 * i + tetFaces[j][k]]];
          p = TMP[0].subVectors(p, center).multiplyScalar(scale).add(center);
          bm.verts.push(new THREE.Vector3(p.x, p.y, p.z));
        }
      }
    }
  }

  let nr = 0;
  for (let i = 0; i < numTets; i++) {
    if (oneFacePerTet) {
      const id0 = faces[4 * i];
      const id1 = faces[4 * i + 1];
      const id2 = faces[4 * i + 2];
      const id3 = faces[4 * i + 3];
      bm.faces.push([id0, id1, id2, id3]);
    }
    else {
      for (let j = 0; j < 4; j++) {
        bm.faces.push([nr, nr + 1, nr + 2]);
        nr = nr + 3;
      }
    }
  }
  return bm;
};