给定三角网格模型是实心的、简单多面体的、密度均匀的。三角面的顶点顺序为逆时针（法向量指向外面）。计算其物理量，如体积，质心和转动惯性张量等。

方法

思想是把模型分解为多个四面体的组合。计算每个四面体的物理量，最后加起来就是整个模型的物理量。 取任意一点为参考点，与三角形网格的顶点组成四面体。以体积为例，因为三角形的顶点顺序为逆时针，所以统一用公式计算时，有的四面体的体积为正，有的为负，最终加起来就是整个模型体积。数学证明就是微积分中的高斯散度定理，用面积分计算体积分。

整个模型的体积分转化为模型表面的每个三角面的面积分之和

一般而言，需要计算体积分：

$V$为模型所占据的空间。$p(x, y, z)$一般为$1, x, y, z, x^2, y^2, z^2, xy, xz, yz$。 由高斯散度定理，有：

$S$为模型的闭合表面。$\mathbf{N}$为指向外面的单位法向量。$\mathbf{F}$为一个向量场，其散度为$p$，即 $\mathbf{\nabla} \cdot \mathbf{F}=p$

模型表面有三角网格表示，其表面积分为所有三角面$\mathcal{F}$的和：

即每个三角面的面积分加起来。

每个三角面的面积分转化为四面体的体积分

可以直接计算每个三角面的面积分，然后求和，就得到需要的物理量。这里将每个三角面转化为四面体，求体积分。证明最开始的思路。

构建四面体

任取一点$\mathbf{O}$作为参考点，与三角面$\mathcal{F}$构成一个四面体${V_\mathcal{F}}$，用高斯散度定理，有三角面的面积分为四面体的体积分减去新构建出来的三个面的面积分。

由于，模型的整个三角面网格构成一个密闭几何网格(watertight mesh)，每一条边对应两个面。所以每一条边总会重复一次，该四面体总有邻接的四面体，新构建出来的三个面的面积分会被邻接的的四面体以相反的法向量再计算一次，求和后就会抵消掉。以正方体的一个面为例，其顶点$a,b,c,d$，组成两个三角面$abc$，$acd$，则其构建的四面体为$abco$，$acdo$，其中面$aco$为共同的面，法向量相反，所以面积分抵消：

所以有：

即每个三角面与参考点构建四面体，计算四面体的体积分，求和。即得到整个模型的体积分。

四面体的体积分计算

对于四面体的体积分，先计算规范化的四面体的体积分，再将其变换到实际的四面体。可以预计算一部分结果，存储起来，方便计算机程序运算。

规范化的四面体由四个点组成${\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3}={(0,0,0), (1,0,0), (0, 1, 0), (0, 0, 1)}$

单位密度均匀分布的刚体的惯性张量为：

有：

考虑其对称性，有：

再将其变换到实际的四面体$\mathbf{C}^{\prime}$，可以分解为旋转放缩$\mathbf{A}$和平移$\Delta \mathbf{x}$两个部分，有：

$\bar{\mathbf{x}}$ 是质心。 对于变换$\mathbf{A}$，有：

$\mathbf{w}_i$为目标四面体，$\mathbf{w}_0$为任意选取的参考点，$\mathbf{w}_1, \mathbf{w}_2, \mathbf{w}_3$为三角面的顶点。 一般，计算以质心为原点的惯性张量，所以：

实现时，参考点$\mathbf{w}_0$直接取质心$\bar{\mathbf{x}}$，省去平移操作。

实现

惯性张量：

void Rigidbody::ComputeInertiaMoment()
{
	// modified from https://github.com/jrouwe/JoltPhysics/Jolt/Physics/Collision/Shape/ConvexHullShape.cpp

	// Calculate covariance matrix
	// See: 
	// - Why the inertia tensor is the inertia tensor - Jonathan Blow (http://number-none.com/blow/inertia/deriving_i.html)
	// - How to find the inertia tensor (or other mass properties) of a 3D solid body represented by a triangle mesh (Draft) - Jonathan Blow, Atman J Binstock (http://number-none.com/blow/inertia/bb_inertia.doc)
	Matrix3f covariance_canonical;
	covariance_canonical << 1.0f / 60.0f, 1.0f / 120.0f, 1.0f / 120.0f,
		1.0f / 120.0f, 1.0f / 60.0f, 1.0f / 120.0f, 
		1.0f / 120.0f, 1.0f / 120.0f, 1.0f / 60.0f;
	Matrix3f covariance_matrix = Matrix3f::Zero(3, 3);

	for (size_t i = 0; i < geometry_.indices.size(); i += 3)
	{
		Vector3f v1 = geometry_.vertices[geometry_.indices[i]].position - centerOfMass_;
		Vector3f v2 = geometry_.vertices[geometry_.indices[i+1]].position - centerOfMass_;
		Vector3f v3 = geometry_.vertices[geometry_.indices[i+2]].position - centerOfMass_;

        // Affine transform that transforms a unit tetrahedon (with vertices (0, 0, 0), (1, 0, 0), (0, 1, 0) and (0, 0, 1) to this tetrahedron
		Matrix3f a;
		a << v1[0], v2[0], v3[0],
        v1[1], v2[1], v3[1],
        v1[2], v2[2], v3[2];

        // Calculate covariance matrix for this tetrahedron
        float det_a = a.determinant();
        Matrix3f c = det_a * (a * covariance_canonical * a.transpose());

        // Add it
        covariance_matrix += c;
	}

	// Calculate inertia matrix assuming density is 1, note that element (3, 3) is garbage
	inertiaMoment = Matrix3f::Identity() * (covariance_matrix(0, 0) + covariance_matrix(1, 1) + covariance_matrix(2, 2)) - covariance_matrix;
	if (abs(volume_) > FLT_MIN && invMass > FLT_MIN)
	{
		inertiaMoment = inertiaMoment / volume_ / invMass;
	}
	else
	{
		inertiaMoment.setZero();
	}
}

质心和体积：

void Rigidbody::ComputeCenterOfMassAndVolume()
{
	// modified from // https://github.com/jrouwe/JoltPhysics/Jolt/Geometry/ConvexHullBuilder.cpp#L1256

	Vector3f v4 = std::accumulate(geometry_.vertices.cbegin(), geometry_.vertices.cend(), Vertex{}, [](const Vertex &a, const Vertex &b)
                               { Vertex vertex; vertex.position = a.position + b.position; return vertex; }).position/geometry_.vertices.size();

	// Calculate mass and center of mass of this convex hull by summing all tetrahedrons
	volume_ = 0.0f;
	centerOfMass_.setZero();
	for (size_t i = 0; i < geometry_.indices.size(); i += 3)
	{
		Vector3f v1 = geometry_.vertices[geometry_.indices[i]].position;
		Vector3f v2 = geometry_.vertices[geometry_.indices[i+1]].position;
		Vector3f v3 = geometry_.vertices[geometry_.indices[i+2]].position;

        // Calculate center of mass and mass of this tetrahedron,
        // see: https://en.wikipedia.org/wiki/Tetrahedron#Volume
        float volume_tetrahedron = (v1 - v4).dot((v2 - v4).cross(v3 - v4)); // Needs to be divided by 6, postpone this until the end of the loop

		auto vt1 = v1 - v4;
		auto vt2 = v2 - v4;
		auto vt3 = v3 - v4;
		auto vt4 = (vt2).cross(vt3);
		auto vt5 = (vt1).dot(vt4);

        Vector3f center_of_mass_tetrahedron = v1 + v2 + v3 + v4; // Needs to be divided by 4, postpone this until the end of the loop

        // Accumulate results
        volume_ += volume_tetrahedron;
        centerOfMass_ += volume_tetrahedron * center_of_mass_tetrahedron;
	}

	// Calculate center of mass, fall back to average point in case there is no volume (everything is on a plane in this case)
	if (volume_ > FLT_EPSILON)
		centerOfMass_ /= 4.0f * volume_;
	else
		centerOfMass_ = v4;

	volume_ /= 6.0f;
}

测试

1. 密度为1，半径为5的球的惯性张量：

• 理论值： $\begin{bmatrix}{\frac {2}{5}}mr^{2}&0&0\\0&{\frac {2}{5}}mr^{2}&0\\0&0&{\frac {2}{5}}mr^{2}\end{bmatrix}$

• 实验：

2. 宽为10，高为5，深为3的长方体的惯性张量：

• 理论值：$\begin{bmatrix}{\frac {1}{12}}m(h^{2}+d^{2})&0&0\\0&{\frac {1}{12}}m(w^{2}+d^{2})&0\\0&0&{\frac {1}{12}}m(w^{2}+h^{2})\end{bmatrix}$

• 实验：

3. TODO: 增加测试样例

参考

1. David Eberly, Polyhedral Mass Properties (Revisited)
2. Jonathan Blow, Atman J Binstock，How to find the inertia tensor (or other mass properties) of a 3D solid body represented by a triangle mesh
3. github.com/jrouwe/Jolt…
4. github.com/davideberly…