0 Introduction
This document explores computational approaches for physics-based simulations in game development. While these methods are commonly used in engineering and scientific computing, their adaptation for real-time interactive applications presents unique challenges.
1 Energy-Based Physics
Physics simulations fundamentally rely on energy conservation principles. Systems possess kinetic energy from motion and potential energy stored in deformations. These energy forms continuously transform in to one another through forces.
Forces emerge as the mechanism driving energy transformations. Consider a stretched spring at rest - it contains potential energy due to its deformed state. The spring exerts forces to return to its equilibrium position, thereby reducing potential energy.
Newton's third law reveals that forces arise from interactions between objects. While the net force across a closed system may be zero, individual objects within the system experience forces that cause acceleration and motion.
In simulation contexts:
- Mass-spring systems model connections between particles using springs defined by energy
- Finite element systems represent deformable objects where energy defines shape
- Finite volume methods simulate volumetric deformations through energy-based constraints
These approaches differ in their geometric representations but share the common principle of defining material behavior through energy functions.
2 Tetrahedral Mesh Representation
2.1 Element Data Format
The element file describes tetrahedral connectivity:
/*
First line: <# of tetrahedra> <nodes per tetrahedron> <# of attributes>
Remaining lines list of # of tetrahedra:
<tetrahedron #> <node> <node> <node> <node> ... [attributes]
*/
1389 4 0
0 1 75 175 52
1 1 75 52 360
0 1 11 46 174
0 1 11 174 57
2.2 Node Data Format
The node file contains vertex positions:
/*
First line: <# of points> <dimension (must be 3)> <# of attributes> <# of boundary markers (0 or 1)>
Remaining lines list # of points:
<point #> <x> <y> <z> [attributes] [boundary marker]
*/
400 3 0 1
1 6.5 1 7.5430812757201648 1
2 4.4456243310934527 5.3728065092415962 8.3129834982130344 0
3 4.1326492294372636 8.1243529245273844 4.2828387696230079 0
4 4.75 7.7548808372563904 7.6928224992291403 0
3 Implementation
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using System;
using System.IO;
public class PhysicsSimulation : MonoBehaviour
{
// Simulation parameters
private float timeStep = 0.003f;
private float particleMass = 1.0f;
private float primaryStiffness = 20000.0f;
private float secondaryStiffness = 5000.0f;
private float dampingFactor = 0.999f;
// Mesh data
private int[] tetrahedra;
private int tetCount;
// Physics state
private Vector3[] forces;
private Vector3[] velocities;
private Vector3[] positions;
private int vertexCount;
// Deformation tracking
private Matrix4x4[] inverseDeformationMatrices;
// Smoothing buffers
private Vector3[] velocitySums;
private int[] neighborCounts;
private MatrixDecomposition matrixSolver = new MatrixDecomposition();
void Initialize()
{
// Load tetrahedral mesh from files
{
string fileContent = File.ReadAllText("Assets/house2.ele");
string[] dataLines = fileContent.Split(new char[] { ' ', '\t', '\r', '\n' }, StringSplitOptions.RemoveEmptyEntries);
tetCount = int.Parse(dataLines[0]);
tetrahedra = new int[tetCount * 4];
for (int tetIndex = 0; tetIndex < tetCount; tetIndex++)
{
tetrahedra[tetIndex * 4] = int.Parse(dataLines[tetIndex * 5 + 4]) - 1;
tetrahedra[tetIndex * 4 + 1] = int.Parse(dataLines[tetIndex * 5 + 5]) - 1;
tetrahedra[tetIndex * 4 + 2] = int.Parse(dataLines[tetIndex * 5 + 6]) - 1;
tetrahedra[tetIndex * 4 + 3] = int.Parse(dataLines[tetIndex * 5 + 7]) - 1;
}
}
{
string fileContent = File.ReadAllText("Assets/house2.node");
string[] dataLines = fileContent.Split(new char[] { ' ', '\t', '\r', '\n' }, StringSplitOptions.RemoveEmptyEntries);
vertexCount = int.Parse(dataLines[0]);
positions = new Vector3[vertexCount];
for (int i = 0; i < vertexCount; i++)
{
positions[i].x = float.Parse(dataLines[i * 5 + 5]) * 0.4f;
positions[i].y = float.Parse(dataLines[i * 5 + 6]) * 0.4f;
positions[i].z = float.Parse(dataLines[i * 5 + 7]) * 0.4f;
}
// Center the model
Vector3 centroid = Vector3.zero;
for (int i = 0; i < vertexCount; i++) centroid += positions[i];
centroid /= vertexCount;
for (int i = 0; i < vertexCount; i++)
{
positions[i] -= centroid;
float temp = positions[i].y;
positions[i].y = positions[i].z;
positions[i].z = temp;
}
}
// Create render mesh
Vector3[] renderVertices = new Vector3[tetCount * 12];
int vertexIndex = 0;
for (int tetIndex = 0; tetIndex < tetCount; tetIndex++)
{
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 2]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 1]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 3]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 2]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 1]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 3]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 1]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 2]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 3]];
}
int[] triangles = new int[tetCount * 12];
for (int tri = 0; tri < tetCount * 4; tri++)
{
triangles[tri * 3] = tri * 3;
triangles[tri * 3 + 1] = tri * 3 + 1;
triangles[tri * 3 + 2] = tri * 3 + 2;
}
Mesh mesh = GetComponent<MeshFilter>().mesh;
mesh.vertices = renderVertices;
mesh.triangles = triangles;
mesh.RecalculateNormals();
// Initialize physics arrays
velocities = new Vector3[vertexCount];
forces = new Vector3[vertexCount];
velocitySums = new Vector3[vertexCount];
neighborCounts = new int[vertexCount];
// Precompute inverse deformation matrices
inverseDeformationMatrices = new Matrix4x4[tetCount];
for (int tetIndex = 0; tetIndex < tetCount; tetIndex++)
{
inverseDeformationMatrices[tetIndex] = ConstructDeformationMatrix(tetIndex).inverse;
}
}
Matrix4x4 ConstructDeformationMatrix(int tetIndex)
{
Matrix4x4 deformationMatrix = Matrix4x4.zero;
Vector3 p1 = positions[tetrahedra[tetIndex * 4 + 1]] - positions[tetrahedra[tetIndex * 4]];
Vector3 p2 = positions[tetrahedra[tetIndex * 4 + 2]] - positions[tetrahedra[tetIndex * 4]];
Vector3 p3 = positions[tetrahedra[tetIndex * 4 + 3]] - positions[tetrahedra[tetIndex * 4]];
deformationMatrix.SetColumn(0, new Vector4(p1.x, p1.y, p1.z, 0.0f));
deformationMatrix.SetColumn(1, new Vector4(p2.x, p2.y, p2.z, 0.0f));
deformationMatrix.SetColumn(2, new Vector4(p3.x, p3.y, p3.z, 0.0f));
deformationMatrix.SetColumn(3, new Vector4(0.0f, 0.0f, 0.0f, 1.0f));
return deformationMatrix;
}
void PhysicsStep()
{
// Apply impulse force
if (Input.GetKeyDown(KeyCode.Space))
{
for (int i = 0; i < vertexCount; i++)
velocities[i].y += 0.2f;
}
// Apply gravity
for (int i = 0; i < vertexCount; i++)
{
forces[i] = new Vector3(0, -9.8f, 0) * particleMass;
}
// Compute elastic forces
for (int tetIndex = 0; tetIndex < tetCount; tetIndex++)
{
// Calculate deformation gradient
Vector3 x1 = positions[tetrahedra[tetIndex * 4 + 1]] - positions[tetrahedra[tetIndex * 4]];
Vector3 x2 = positions[tetrahedra[tetIndex * 4 + 2]] - positions[tetrahedra[tetIndex * 4]];
Vector3 x3 = positions[tetrahedra[tetIndex * 4 + 3]] - positions[tetrahedra[tetIndex * 4]];
Matrix4x4 deformationMatrix = new Matrix4x4();
deformationMatrix.SetColumn(0, new Vector4(x1.x, x1.y, x1.z, 0.0f));
deformationMatrix.SetColumn(1, new Vector4(x2.x, x2.y, x2.z, 0.0f));
deformationMatrix.SetColumn(2, new Vector4(x3.x, x3.y, x3.z, 0.0f));
deformationMatrix.SetColumn(3, new Vector4(0.0f, 0.0f, 0.0f, 1.0f));
Matrix4x4 deformationGradient = deformationMatrix * inverseDeformationMatrices[tetIndex];
// Calculate Green-Lagrange strain
Matrix4x4 strainTensor = ScaleMatrix(0.5f, SubtractMatrices(
MultiplyMatrices(deformationGradient.transpose, deformationGradient),
Matrix4x4.identity));
// Compute second Piola-Kirchhoff stress
Matrix4x4 stressTensor = AddMatrices(
ScaleMatrix(2 * secondaryStiffness, strainTensor),
ScaleMatrix(primaryStiffness * MatrixTrace(strainTensor), Matrix4x4.identity));
Matrix4x4 firstPiolaKirchhoff = MultiplyMatrices(deformationGradient, stressTensor);
// Calculate elastic forces
Matrix4x4 elasticForces = ScaleMatrix(-1.0f / (6.0f * inverseDeformationMatrices[tetIndex].determinant),
MultiplyMatrices(firstPiolaKirchhoff, inverseDeformationMatrices[tetIndex].transpose));
// Distribute forces to vertices
Vector3 f1 = new Vector3(elasticForces.GetColumn(0).x, elasticForces.GetColumn(0).y, elasticForces.GetColumn(0).z);
Vector3 f2 = new Vector3(elasticForces.GetColumn(1).x, elasticForces.GetColumn(1).y, elasticForces.GetColumn(1).z);
Vector3 f3 = new Vector3(elasticForces.GetColumn(2).x, elasticForces.GetColumn(2).y, elasticForces.GetColumn(2).z);
Vector3 f0 = -(f1 + f2 + f3);
forces[tetrahedra[tetIndex * 4]] += f0;
forces[tetrahedra[tetIndex * 4 + 1]] += f1;
forces[tetrahedra[tetIndex * 4 + 2]] += f2;
forces[tetrahedra[tetIndex * 4 + 3]] += f3;
}
// Update positions and velocities
for (int i = 0; i < vertexCount; i++)
{
// Semi-implicit Euler integration
velocities[i] += forces[i] * timeStep / particleMass;
velocities[i] *= dampingFactor;
positions[i] += velocities[i] * timeStep;
// Floor collision
if (positions[i].y < -3)
{
positions[i].y = -3;
velocities[i] *= 0.5f;
}
}
}
void Update()
{
// Multiple physics substeps for stability
for (int step = 0; step < 10; step++)
PhysicsStep();
// Update render mesh
Vector3[] renderVertices = new Vector3[tetCount * 12];
int vertexIndex = 0;
for (int tetIndex = 0; tetIndex < tetCount; tetIndex++)
{
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 2]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 1]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 3]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 2]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 1]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 3]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 1]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 2]];
renderVertices[vertexIndex++] = positions[tetrahedra[tetIndex * 4 + 3]];
}
Mesh mesh = GetComponent<MeshFilter>().mesh;
mesh.vertices = renderVertices;
mesh.RecalculateNormals();
}
// Matrix utility methods
private Matrix4x4 SubtractMatrices(Matrix4x4 left, Matrix4x4 right)
{
Matrix4x4 result = Matrix4x4.zero;
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++)
result[i, j] = left[i, j] - right[i, j];
result[3, 3] = 1;
return result;
}
private Matrix4x4 AddMatrices(Matrix4x4 left, Matrix4x4 right)
{
Matrix4x4 result = Matrix4x4.zero;
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++)
result[i, j] = left[i, j] + right[i, j];
result[3, 3] = 1;
return result;
}
private Matrix4x4 ScaleMatrix(float scalar, Matrix4x4 matrix)
{
Matrix4x4 result = Matrix4x4.zero;
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++)
result[i, j] = scalar * matrix[i, j];
result[3, 3] = 1;
return result;
}
private float MatrixTrace(Matrix4x4 matrix)
{
return matrix[0, 0] + matrix[1, 1] + matrix[2, 2];
}
private Matrix4x4 MultiplyMatrices(Matrix4x4 left, Matrix4x4 right)
{
Matrix4x4 result = Matrix4x4.zero;
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++)
for (int k = 0; k < 4; k++)
result[i, j] += left[i, k] * right[k, j];
return result;
}
}
// Helper class for matrix operations (placeholder)
public class MatrixDecomposition
{
// Implementation would include SVD and other decompositions
}