Sim/Part.cs
using System;
using System.Collections.Generic;

namespace SquishyTerry.Sim;

/// <summary>
/// Small soft body glued to the face: an eye, the nose or an ear. Built from a clipped level-3
/// icosphere and anchored to one face vertex plus that feature patch's shape-matching rotation.
/// Original: 172-byte class, ctor f_hc, update f_gc.
/// </summary>
public sealed class Part
{
	public Vec3[] Pos, PrevPos, Pred, Vel, PrevVel, Normals;
	public int[] Indices;
	/// <summary>Indices plus a fan cap over the base rim, used for volume.</summary>
	public int[] ClosedIndices;
	public float[] InvMass;
	/// <summary>Normalized height 0..1 above the cut plane (before scaling).</summary>
	public float[] Height;
	public Vec3[] RestLocal;
	/// <summary>Optional per-vertex shading mask (e.g. nostrils), or null.</summary>
	public float[] Mask;
	/// <summary>Rest shape before it was rotated onto the anchor normal (+Z out of the face, +Y up).</summary>
	public Vec3[] RestUnaligned;
	/// <summary>Rotation taking <see cref="RestUnaligned"/> to <see cref="RestLocal"/>.</summary>
	public Quat Align;
	readonly DistanceConstraints _distance;
	// Per-vertex soft blend weight f (constant, from Height) and 1 - f; goal-position scratch.
	readonly float[] _blend, _blendKeep;
	readonly Vec3[] _goal;

	public int VertexCount => Pos.Length;

	/// <param name="shape">Optional sculpt applied to the unit sphere before clipping (null in the original).</param>
	/// <param name="detail">Icosphere subdivisions (3 in the original).</param>
	/// <param name="mask">Optional per-vertex shading mask, evaluated on the sculpted unit sphere (see <see cref="Mask"/>).</param>
	public Part( Vec3 anchor, Vec3 normal, Vec3 scale, float cutZ, float shear, float rotZ, Func<Vec3, Vec3> shape = null, int detail = 3, Func<Vec3, float> mask = null )
	{
		var pos = new List<Vec3>();
		var idx = new List<int>();
		Geometry.Icosphere( detail, pos, idx );
		if ( shape is not null )
		{
			for ( int i = 0; i < pos.Count; i++ )
				pos[i] = shape( pos[i] );
		}

		// --- Clip the sphere at z >= cutZ (f_cd / f_bd / f_mc) ---
		var clipped = new List<bool>( new bool[pos.Count] );
		var cutVerts = new Dictionary<long, int>();
		var outIdx = new List<int>();

		int CutVertex( int v1, int v2, float d1, float d2 )
		{
			int lo = Math.Min( v1, v2 ), hi = Math.Max( v1, v2 );
			long key = ((long)lo << 32) | (uint)hi;
			if ( cutVerts.TryGetValue( key, out int existing ) ) return existing;
			float dLo = MathF.Abs( v1 > v2 ? d2 : d1 );
			float dHi = MathF.Abs( v1 > v2 ? d1 : d2 );
			float r = dLo / (dLo + dHi);
			pos.Add( pos[hi] * r + pos[lo] * (1.0f - r) );
			clipped.Add( true );
			cutVerts[key] = pos.Count - 1;
			return pos.Count - 1;
		}

		void OneIn( int a, int b, int c, float da, float db, float dc )
		{
			if ( !(da > 0 && db < 0 && dc < 0) ) return;
			int k = CutVertex( a, b, da, db );
			int h = CutVertex( a, c, da, dc );
			outIdx.Add( a ); outIdx.Add( k ); outIdx.Add( h );
		}

		void TwoIn( int a, int b, int c, float da, float db, float dc )
		{
			if ( !(da > 0 && db > 0 && dc < 0) ) return;
			int n = CutVertex( a, c, da, dc );
			int j = CutVertex( b, c, db, dc );
			outIdx.Add( a ); outIdx.Add( b ); outIdx.Add( j );
			outIdx.Add( a ); outIdx.Add( j ); outIdx.Add( n );
		}

		for ( int t = 0; t + 2 < idx.Count; t += 3 )
		{
			int a = idx[t], b = idx[t + 1], c = idx[t + 2];
			float da = pos[a].Z - cutZ, db = pos[b].Z - cutZ, dc = pos[c].Z - cutZ;
			if ( da >= 0 && db >= 0 && dc >= 0 )
			{
				outIdx.Add( a ); outIdx.Add( b ); outIdx.Add( c );
				continue;
			}
			OneIn( a, b, c, da, db, dc );
			OneIn( b, c, a, db, dc, da );
			OneIn( c, a, b, dc, da, db );
			TwoIn( a, b, c, da, db, dc );
			TwoIn( b, c, a, db, dc, da );
			TwoIn( c, a, b, dc, da, db );
		}

		// --- Drop unused vertices, keeping order (f_pq) ---
		var used = new bool[pos.Count];
		foreach ( int i in outIdx ) used[i] = true;
		var remap = new int[pos.Count];
		var newPos = new List<Vec3>();
		var newClipped = new List<bool>();
		for ( int i = 0; i < pos.Count; i++ )
		{
			if ( !used[i] ) continue;
			remap[i] = newPos.Count;
			newPos.Add( pos[i] );
			newClipped.Add( clipped[i] );
		}
		Indices = new int[outIdx.Count];
		for ( int i = 0; i < outIdx.Count; i++ ) Indices[i] = remap[outIdx[i]];
		int n = newPos.Count;
		var p = newPos.ToArray();

		if ( mask is not null )
		{
			Mask = new float[n];
			for ( int i = 0; i < n; i++ ) Mask[i] = mask( p[i] );
		}

		// --- Shape it ---
		for ( int i = 0; i < n; i++ )
		{
			float z = p[i].Z - cutZ;
			p[i] = new Vec3( p[i].X + z * shear, p[i].Y, z );
		}

		Height = new float[n];
		float hs = 1.0f / (1.0f - cutZ);
		for ( int i = 0; i < n; i++ ) Height[i] = p[i].Z * hs;

		float sinZ = MathF.Sin( rotZ * 0.5f ), cosZ = MathF.Cos( rotZ * 0.5f );
		float two = sinZ * 2.0f;
		for ( int i = 0; i < n; i++ )
		{
			float x = p[i].X * scale.X, y = p[i].Y * scale.Y, z = p[i].Z * scale.Z;
			p[i] = new Vec3( x - two * (y * cosZ + x * sinZ), two * (x * cosZ - y * sinZ) + y, z );
		}

		// --- Base rim loop and closed index list (f_ry) ---
		var next = new Dictionary<int, int>();
		for ( int t = 0; t + 2 < Indices.Length; t += 3 )
		{
			int a = Indices[t], b = Indices[t + 1], c = Indices[t + 2];
			if ( newClipped[a] && newClipped[b] ) next[a] = b;
			if ( newClipped[b] && newClipped[c] ) next[b] = c;
			if ( newClipped[c] && newClipped[a] ) next[c] = a;
		}
		var rim = new List<int>();
		int start = newClipped.IndexOf( true );
		rim.Add( start );
		while ( true )
		{
			next.TryGetValue( rim[^1], out int nx );
			if ( nx == rim[0] || rim.Count > n ) break;
			rim.Add( nx );
		}
		var closed = new List<int>( Indices );
		for ( int j = 1; j < rim.Count - 1; j++ )
		{
			closed.Add( rim[0] ); closed.Add( rim[j] ); closed.Add( rim[j + 1] );
		}
		ClosedIndices = closed.ToArray();

		InvMass = new float[n];
		for ( int i = 0; i < n; i++ ) InvMass[i] = p[i].Z <= 0.02f ? 0.0f : 1.0f;

		// --- Rotate +Z onto the anchor normal ---
		float xy = MathF.Sqrt( normal.X * normal.X + normal.Y * normal.Y );
		float half = MathF.Asin( xy ) * 0.5f;
		float k = 1.0f / xy * MathF.Sin( half );
		// The original divides by |n.xy|; guard the straight-ahead case it never hit.
		var align = xy > 1e-6f ? new Quat( -normal.Y * k, normal.X * k, 0, MathF.Cos( half ) ) : Quat.Identity;
		Align = align;
		RestUnaligned = Geometry.Copy( p );
		for ( int i = 0; i < n; i++ ) p[i] = align.Rotate( p[i] );

		RestLocal = Geometry.Copy( p );
		Pos = new Vec3[n];
		for ( int i = 0; i < n; i++ ) Pos[i] = p[i] + anchor;
		PrevPos = Geometry.Copy( Pos );
		Pred = new Vec3[n];
		Vel = new Vec3[n];
		PrevVel = new Vec3[n];
		Normals = new Vec3[n];
		Geometry.ComputeNormals( Pos, Indices, Normals );

		_distance = new DistanceConstraints( Geometry.BuildEdges( Pos, Indices, 0.85f ) );
		_goal = new Vec3[n];
		_blend = new float[n];
		_blendKeep = new float[n];
		for ( int i = 0; i < n; i++ )
		{
			float s = Math.Clamp( Height[i] * 2.0f, 0.0f, 1.0f );
			float f = s * s * 0.7f * (3.0f - s * 2.0f);
			_blend[i] = f;
			_blendKeep[i] = 1.0f - f;
		}
		// The original also builds bending constraints here but never solves them.
	}

	/// <param name="anchor">Face vertex position the part sits on.</param>
	/// <param name="normal">Face vertex normal.</param>
	/// <param name="rotation">Rotation of the underlying shape-matching patch.</param>
	/// <param name="anchorVelocity">Face vertex velocity; above 0.05 the part goes nearly rigid.</param>
	public void Update( Vec3 anchor, Vec3 normal, Quat rotation, Vec3 anchorVelocity )
	{
		// Vec3 math written out by hand (the editor build is unoptimized), in the original evaluation
		// order. Loops that only touch vertex i are fused; goal = anchor + rotation * RestLocal is
		// computed once per vertex.
		Vec3[] pos = Pos, prevPos = PrevPos, pred = Pred, vel = Vel, prevVel = PrevVel, restLocal = RestLocal, goal = _goal;
		float[] invMass = InvMass;
		int n = pos.Length;
		float ax = anchor.X, ay = anchor.Y, az = anchor.Z;
		float qx = rotation.X, qy = rotation.Y, qz = rotation.Z, qw = rotation.W;
		for ( int i = 0; i < n; i++ )
		{
			ref Vec3 vl = ref vel[i];
			vl.X *= 0.98f; vl.Y *= 0.98f; vl.Z *= 0.98f;

			ref Vec3 r = ref restLocal[i];
			float tx = 2f * (qy * r.Z - qz * r.Y);
			float ty = 2f * (qz * r.X - qx * r.Z);
			float tz = 2f * (qx * r.Y - qy * r.X);
			ref Vec3 g = ref goal[i];
			g.X = ax + (r.X + qw * tx + (qy * tz - qz * ty));
			g.Y = ay + (r.Y + qw * ty + (qz * tx - qx * tz));
			g.Z = az + (r.Z + qw * tz + (qx * ty - qy * tx));

			ref Vec3 pr = ref pred[i];
			if ( invMass[i] == 0.0f )
			{
				pr = g;
			}
			else
			{
				ref Vec3 p = ref pos[i];
				pr.X = p.X + vl.X; pr.Y = p.Y + vl.Y; pr.Z = p.Z + vl.Z;
			}
		}

		Geometry.SolveDistance( _distance, pred, invMass );

		bool rigid = MathF.Sqrt( anchorVelocity.X * anchorVelocity.X + anchorVelocity.Y * anchorVelocity.Y + anchorVelocity.Z * anchorVelocity.Z ) > 0.05f;
		float[] blend = _blend, blendKeep = _blendKeep;
		float nx = normal.X, ny = normal.Y, nz = normal.Z;
		for ( int i = 0; i < n; i++ )
		{
			ref Vec3 g = ref goal[i];
			ref Vec3 pr = ref pred[i];
			if ( rigid )
			{
				pr.X = g.X * 0.9f + pr.X * 0.1f;
				pr.Y = g.Y * 0.9f + pr.Y * 0.1f;
				pr.Z = g.Z * 0.9f + pr.Z * 0.1f;
			}
			else
			{
				float f = blend[i], keep = blendKeep[i];
				pr.X = g.X * keep + pr.X * f;
				pr.Y = g.Y * keep + pr.Y * f;
				pr.Z = g.Z * keep + pr.Z * f;
			}

			// Keep everything above the face's tangent plane.
			float d = (pr.X - ax) * nx + (pr.Y - ay) * ny + (pr.Z - az) * nz;
			if ( d < 0.0f && invMass[i] != 0.0f )
			{
				pr.X -= nx * d; pr.Y -= ny * d; pr.Z -= nz * d;
			}

			ref Vec3 p = ref pos[i];
			ref Vec3 vl = ref vel[i];
			prevVel[i] = vl;
			vl.X = pr.X - p.X; vl.Y = pr.Y - p.Y; vl.Z = pr.Z - p.Z;
			prevPos[i] = p;
			p = pr;
		}
		Geometry.ComputeNormals( pos, Indices, Normals );
	}

	public float ComputeVolume() => Metrics.Volume( Pos, ClosedIndices );

	public float ComputeAverageSpeed() => Metrics.AverageSpeed( Vel );
}