Code/HumanoidRetargeter/Core/Formats/Mmd/Bullet/BulletDynamics/JacobianEntry.cs

A C# port of a Bullet physics Jacobian entry used by constraint solvers. It computes linear/angular Jacobian terms, diagonal and off-diagonal entries, and relative velocity for constraints between one or two rigid bodies.

Native Interop
/*
 * C# / XNA  port of Bullet (c) 2011 Mark Neale <[email protected]>
 *
 * Bullet Continuous Collision Detection and Physics Library
 * Copyright (c) 2003-2008 Erwin Coumans  http://www.bulletphysics.com/
 *
 * This software is provided 'as-is', without any express or implied warranty.
 * In no event will the authors be held liable for any damages arising from
 * the use of this software.
 *
 * Permission is granted to anyone to use this software for any purpose,
 * including commercial applications, and to alter it and redistribute it
 * freely, subject to the following restrictions:
 *
 * 1. The origin of this software must not be misrepresented; you must not
 *    claim that you wrote the original software. If you use this software
 *    in a product, an acknowledgment in the product documentation would be
 *    appreciated but is not required.
 * 2. Altered source versions must be plainly marked as such, and must not be
 *    misrepresented as being the original software.
 * 3. This notice may not be removed or altered from any source distribution.
 */
// Adapted from BulletXNA 73b2e1b9: namespace isolation, no XNA rendering, one type per file.
#nullable disable
#pragma warning disable CS1591, CS0168, CS0219, CS0414, CS0660, CS0661, CS0659, CS0693, CS8981, CS0108, CS0114
using System.Diagnostics;
using HumanoidRetargeter.Core.Formats.Mmd.Bullet.LinearMath;

namespace HumanoidRetargeter.Core.Formats.Mmd.Bullet.BulletDynamics;

internal class JacobianEntry
{
	public JacobianEntry()
	{
	}

	//constraint between two different rigidbodies
	public JacobianEntry(IndexedBasisMatrix world2A, IndexedBasisMatrix world2B, IndexedVector3 rel_pos1, IndexedVector3 rel_pos2, IndexedVector3 jointAxis, IndexedVector3 inertiaInvA, float massInvA, IndexedVector3 inertiaInvB, float massInvB) : this(ref world2A, ref world2B, ref rel_pos1, ref rel_pos2, ref jointAxis, ref inertiaInvA, massInvA, ref inertiaInvB, massInvB)
	{
	}

	public JacobianEntry(ref IndexedBasisMatrix world2A, ref IndexedBasisMatrix world2B, ref IndexedVector3 rel_pos1, ref IndexedVector3 rel_pos2, ref IndexedVector3 jointAxis, ref IndexedVector3 inertiaInvA, float massInvA, ref IndexedVector3 inertiaInvB, float massInvB)
	{
		m_linearJointAxis = jointAxis;
		m_aJ = world2A * (rel_pos1.Cross(ref m_linearJointAxis));
		m_bJ = world2B * (rel_pos2.Cross(-m_linearJointAxis));
		m_0MinvJt = inertiaInvA * m_aJ;
		m_1MinvJt = inertiaInvB * m_bJ;
		m_Adiag = massInvA + IndexedVector3.Dot(m_0MinvJt, m_aJ) + massInvB + IndexedVector3.Dot(m_1MinvJt, m_bJ);
		{
		}
	}

	//angular constraint between two different rigidbodies
	public JacobianEntry(IndexedVector3 jointAxis, IndexedBasisMatrix world2A, IndexedBasisMatrix world2B, IndexedVector3 inertiaInvA, IndexedVector3 inertiaInvB) : this(ref jointAxis, ref world2A, ref world2B, ref inertiaInvA, ref inertiaInvB)
	{
	}

	public JacobianEntry(ref IndexedVector3 jointAxis, ref IndexedBasisMatrix world2A, ref IndexedBasisMatrix world2B, ref IndexedVector3 inertiaInvA, ref IndexedVector3 inertiaInvB)
	{
		m_linearJointAxis = IndexedVector3.Zero;
		m_aJ = world2A * jointAxis;
		m_bJ = world2B * -jointAxis;
		m_0MinvJt = inertiaInvA * m_aJ;
		m_1MinvJt = inertiaInvB * m_bJ;
		m_Adiag = IndexedVector3.Dot(m_0MinvJt, m_aJ) + IndexedVector3.Dot(m_1MinvJt, m_bJ);
		{
		}
	}

	public JacobianEntry(IndexedVector3 axisInA, IndexedVector3 axisInB, IndexedVector3 inertiaInvA, IndexedVector3 inertiaInvB) : this(ref axisInA, ref axisInB, ref inertiaInvA, ref inertiaInvB)
	{
	}

	//angular constraint between two different rigidbodies
	public JacobianEntry(ref IndexedVector3 axisInA, ref IndexedVector3 axisInB, ref IndexedVector3 inertiaInvA, ref IndexedVector3 inertiaInvB)
	{
		m_linearJointAxis = IndexedVector3.Zero;
		m_aJ = axisInA;
		m_bJ = -axisInB;
		m_0MinvJt = inertiaInvA * m_aJ;
		m_1MinvJt = inertiaInvB * m_bJ;
		m_Adiag = IndexedVector3.Dot(m_0MinvJt, m_aJ) + IndexedVector3.Dot(m_1MinvJt, m_bJ);
		{
		}
	}

	//constraint on one rigidbody
	public JacobianEntry(IndexedBasisMatrix world2A, IndexedVector3 rel_pos1, IndexedVector3 rel_pos2, IndexedVector3 jointAxis, IndexedVector3 inertiaInvA, float massInvA) : this(ref world2A, ref rel_pos1, ref rel_pos2, ref jointAxis, ref inertiaInvA, massInvA)
	{
	}

	public JacobianEntry(ref IndexedBasisMatrix world2A, ref IndexedVector3 rel_pos1, ref IndexedVector3 rel_pos2, ref IndexedVector3 jointAxis, ref IndexedVector3 inertiaInvA, float massInvA)
	{
		m_linearJointAxis = jointAxis;
		m_aJ = world2A * (rel_pos1.Cross(ref jointAxis));
		m_bJ = world2A * (rel_pos2.Cross(-jointAxis));
		m_0MinvJt = inertiaInvA * m_aJ;
		m_1MinvJt = IndexedVector3.Zero;
		m_Adiag = massInvA + IndexedVector3.Dot(m_0MinvJt, m_aJ);
		{
		}
	}

	public float GetDiagonal()
	{
		return m_Adiag;
	}

	// for two constraints on the same rigidbody (for example vehicle friction)
	public float GetNonDiagonal(JacobianEntry jacB, float massInvA)
	{
		JacobianEntry jacA = this;
		float lin = massInvA * IndexedVector3.Dot(jacA.m_linearJointAxis, jacB.m_linearJointAxis);
		float ang = IndexedVector3.Dot(jacA.m_0MinvJt, jacB.m_aJ);
		return lin + ang;
	}

	// for two constraints on sharing two same rigidbodies (for example two contact points between two rigidbodies)
	public float GetNonDiagonal(JacobianEntry jacB, float massInvA, float massInvB)
	{
		JacobianEntry jacA = this;
		IndexedVector3 lin = jacA.m_linearJointAxis * jacB.m_linearJointAxis;
		IndexedVector3 ang0 = jacA.m_0MinvJt * jacB.m_aJ;
		IndexedVector3 ang1 = jacA.m_1MinvJt * jacB.m_bJ;
		IndexedVector3 lin0 = massInvA * lin;
		IndexedVector3 lin1 = massInvB * lin;
		IndexedVector3 sum = ang0 + ang1 + lin0 + lin1;
		return sum.X + sum.Y + sum.Z;
	}

	public float GetRelativeVelocity(IndexedVector3 linvelA, IndexedVector3 angvelA, IndexedVector3 linvelB, IndexedVector3 angvelB)
	{
		return GetRelativeVelocity(ref linvelA, ref angvelA, ref linvelB, ref angvelB);
	}

	public float GetRelativeVelocity(ref IndexedVector3 linvelA, ref IndexedVector3 angvelA, ref IndexedVector3 linvelB, ref IndexedVector3 angvelB)
	{
		IndexedVector3 linrel = linvelA - linvelB;
		IndexedVector3 angvela = angvelA * m_aJ;
		IndexedVector3 angvelb = angvelB * m_bJ;
		linrel *= m_linearJointAxis;
		angvela += angvelb;
		angvela += linrel;
		float rel_vel2 = angvela.X + angvela.Y + angvela.Z;
		return rel_vel2 + MathUtil.SIMD_EPSILON;
	}

	public IndexedVector3 m_linearJointAxis;
	public IndexedVector3 m_aJ;
	public IndexedVector3 m_bJ;
	public IndexedVector3 m_0MinvJt;
	public IndexedVector3 m_1MinvJt;
	//Optimization: can be stored in the w/last component of one of the vectors
	public float m_Adiag;
}