Use sphere mesh. Add soft body test.

This commit is contained in:
Irlan
2019-04-18 19:44:12 -03:00
parent 695514989e
commit 399a6efc72
10 changed files with 588 additions and 653 deletions

View File

@ -1,158 +0,0 @@
/*
* Copyright (c) 2016-2019 Irlan Robson https://irlanrobson.github.io
*
* 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.
*/
#ifndef MASS_SPRING_H
#define MASS_SPRING_H
class MassSpring : public Test
{
public:
MassSpring()
{
m_x.Set(0.0f, 5.0f, 0.0f);
m_v.SetZero();
m_k = 100.0f;
m_iterations = 0;
}
void Solve(float32 h)
{
// ODE
// f(Y) = dY / dt = [v]
// [-k * x]
// 1. Apply Implicit Euler
// Y(t + h) = Y(t) + h * f( Y(t + h) )
// G( Y(t + h) ) = Y(t + h) - Y(t) - h * f( Y(t + h) ) = 0
// 2. Solve G = 0
// Newton-Raphson Iteration
//
// Y(t + h) =
// Y(t + h)_0 - G( Y(t + h)_0 ) / G'( Y(t + h)_0 ) =
// Y(t + h)_0 - G'( Y(t + h)_0 )^-1 * ( Y(t + h)_0 - Y(t) - h * f( Y(t + h)_0 )
// G'( Y ) = I - h * del_f / del_Y
// del_f / del_Y = [del_f1 / del_x del_f1 / del_v] = [0 I]
// [del_f2 / del_x del_f2 / del_v] [-k * I 0]
// G'( Y ) = [I 0] - [0 h * I] = [I -h * I]
// [0 I] [-h * k * I 0] [h * k * I I]
// Compute Jacobian
b3Mat33 I = b3Mat33_identity;
b3Mat33 A, B, C, D;
A = I;
B = -h * I;
C = h * m_k * I;
D = I;
// Invert
// Block matrix inversion
b3Mat33 invD = b3Inverse(D);
b3Mat33 B_invD = B * invD;
b3Mat33 invJ_A = b3Inverse(A - B_invD * C);
b3Mat33 invJ_B = -invJ_A * B_invD;
b3Mat33 invJ_C = -invD * C * invJ_A;
b3Mat33 invJ_D = invD + invD * C * invJ_A * B_invD;
// Initial guess
b3Vec3 f1 = m_v;
b3Vec3 f2 = -m_k * m_x;
b3Vec3 Y1 = m_x + h * f1;
b3Vec3 Y2 = m_v + h * f2;
const float32 kTol = 0.05f;
const u32 kMaxIterations = 20;
float32 eps0 = 0.0f;
float32 eps1 = B3_MAX_FLOAT;
m_iterations = 0;
while (m_iterations < kMaxIterations && eps1 > kTol * kTol * eps0)
{
// Evaluate f(Y_n-1)
f1 = Y2;
f2 = -m_k * Y1;
// Residual vector
b3Vec3 G1 = Y1 - m_x - h * f1;
b3Vec3 G2 = Y2 - m_v - h * f2;
eps1 = b3Dot(G1, G1) + b3Dot(G2, G2);
// Solve Ax = b
b3Vec3 x1 = invJ_A * G1 + invJ_B * G2;
b3Vec3 x2 = invJ_C * G1 + invJ_D * G2;
Y1 -= x1;
Y2 -= x2;
++m_iterations;
}
// Update state
m_x = Y1;
m_v = Y2;
}
void Step()
{
float32 h = g_testSettings->inv_hertz;
Solve(h);
g_draw->DrawSolidSphere(m_x, 0.25f, b3Color_white);
g_draw->DrawSegment(b3Vec3_zero, m_x, b3Color_white);
g_draw->DrawString(b3Color_white, "Iterations = %u", m_iterations);
float32 E = 0.5f * b3Dot(m_v, m_v);
g_draw->DrawString(b3Color_white, "E = %f", E);
}
static Test* Create()
{
return new MassSpring();
}
// State
b3Vec3 m_x, m_v;
// Stiffness
float32 m_k;
//
u32 m_iterations;
};
#endif

View File

@ -1,142 +0,0 @@
/*
* Copyright (c) 2016-2019 Irlan Robson https://irlanrobson.github.io
*
* 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.
*/
#ifndef SHIRT_H
#define SHIRT_H
class Shirt : public Test
{
public:
Shirt()
{
// Generate 2D mesh
m_shirtGarmentMesh.Set(&m_shirtGarment, 0.1f);
// Create 3D mesh
m_shirtClothMesh.Set(&m_shirtGarmentMesh);
// Create cloth
b3ClothDef def;
def.mesh = &m_shirtClothMesh;
def.density = 0.2f;
def.structural = 1000.0f;
m_cloth = new b3Cloth(def);
// Perform fitting
for (u32 i = 0; i < 3; ++i)
{
b3ClothMeshMesh* front = m_shirtClothMesh.meshes + i;
for (u32 j = 0; j < front->vertexCount; ++j)
{
u32 v = front->startVertex + j;
b3Particle* p = m_cloth->GetVertexParticle(v);
b3Vec3 x = p->GetPosition();
x.z = -1.0f;
p->SetPosition(x);
}
}
for (u32 i = 3; i < 6; ++i)
{
b3ClothMeshMesh* back = m_shirtClothMesh.meshes + i;
for (u32 j = 0; j < back->vertexCount; ++j)
{
u32 v = back->startVertex + j;
b3Particle* p = m_cloth->GetVertexParticle(v);
b3Vec3 x = p->GetPosition();
x.z = 1.0f;
p->SetPosition(x);
}
}
m_cloth->SetGravity(b3Vec3(0.0f, -9.8f, 0.0f));
m_cloth->SetWorld(&m_world);
m_clothDragger = new b3ClothDragger(&m_ray, m_cloth);
}
~Shirt()
{
delete m_clothDragger;
delete m_cloth;
}
void Step()
{
Test::Step();
m_cloth->Step(g_testSettings->inv_hertz);
m_cloth->Draw();
extern u32 b3_clothSolverIterations;
g_draw->DrawString(b3Color_white, "Iterations = %d", b3_clothSolverIterations);
float32 E = m_cloth->GetEnergy();
g_draw->DrawString(b3Color_white, "E = %f", E);
}
void MouseMove(const b3Ray3& pw)
{
Test::MouseMove(pw);
if (m_clothDragger->IsDragging() == true)
{
m_clothDragger->Drag();
}
}
void MouseLeftDown(const b3Ray3& pw)
{
Test::MouseLeftDown(pw);
if (m_clothDragger->IsDragging() == false)
{
m_clothDragger->StartDragging();
}
}
void MouseLeftUp(const b3Ray3& pw)
{
Test::MouseLeftUp(pw);
if (m_clothDragger->IsDragging() == true)
{
m_clothDragger->StopDragging();
}
}
static Test* Create()
{
return new Shirt();
}
b3ShirtGarment m_shirtGarment;
b3GarmentMesh m_shirtGarmentMesh;
b3GarmentClothMesh m_shirtClothMesh;
b3Cloth* m_cloth;
b3ClothDragger* m_clothDragger;
};
#endif

View File

@ -1,99 +0,0 @@
/*
* Copyright (c) 2016-2019 Irlan Robson https://irlanrobson.github.io
*
* 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.
*/
#ifndef PENDULUM_H
#define PENDULUM_H
class SinglePendulum : public Test
{
public:
SinglePendulum()
{
m_g = -10.0f;
m_r = 10.0f;
m_m = 1.0f;
m_I = m_m * m_r * m_r;
// Initial state
m_theta = -0.5f * B3_PI;
m_omega = 0.0f;
}
void Step()
{
float32 h = g_testSettings->inv_hertz;
// Solution (acceleration)
float32 omega_dot = -m_g / m_r * sin(m_theta);
// Integrate acceleration
m_omega += h * omega_dot;
// Integrate velocity
m_theta += h * m_omega;
// Convert from polar coordinates (r, theta) to Cartesian coordinates (x, y)
b3Vec3 c;
c.x = m_r * sin(m_theta);
c.y = m_r * cos(m_theta);
c.z = 0.0f;
g_draw->DrawSolidSphere(c, 1.0f, b3Color_white);
b3Vec3 pole;
pole.SetZero();
g_draw->DrawSegment(pole, c, b3Color_white);
// Kinetic energy
float32 T = 0.5f * m_I * m_omega * m_omega;
// Potential energy
float32 V = -m_m * m_g * m_r * cos(m_theta);
// Lagrangian
float32 L = T - V;
//
g_draw->DrawString(b3Color_white, "T = %f \nV = %f \nL = %f", T, V, L);
}
static Test* Create()
{
return new SinglePendulum();
}
// Gravity
float32 m_g;
// Mass, inertia
float32 m_m, m_I;
// Radial coordinate
float32 m_r;
// The allowable generalized coordinate in polar coordinate frame.
// Only motions satisfying the constraints can be described
// in this frame. Therefore, all solutions satisfy the constraints.
// This is the so called reduced coordinates approach.
float32 m_theta;
// Velocity
float32 m_omega;
};
#endif

View File

@ -0,0 +1,278 @@
/*
* Copyright (c) 2016-2019 Irlan Robson https://irlanrobson.github.io
*
* 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.
*/
#ifndef SOFT_BODY_H
#define SOFT_BODY_H
#include <testbed/framework/sphere_mesh.h>
struct b3SphereClothMesh : public b3ClothMesh
{
b3StackArray<b3Vec3, 256> sphereVertices;
b3StackArray<b3ClothMeshTriangle, 256> sphereTriangles;
b3ClothMeshMesh sphereMesh;
b3SphereClothMesh()
{
smMesh mesh;
smCreateMesh(mesh, 1);
sphereVertices.Resize(mesh.vertices.Count());
for (u32 i = 0; i < mesh.vertices.Count(); ++i)
{
sphereVertices[i] = mesh.vertices[i];
}
sphereTriangles.Resize(mesh.triangleIndices.Count() / 3);
for (u32 i = 0; i < mesh.triangleIndices.Count() / 3; ++i)
{
sphereTriangles[i].v1 = mesh.triangleIndices[3 * i + 0];
sphereTriangles[i].v2 = mesh.triangleIndices[3 * i + 1];
sphereTriangles[i].v3 = mesh.triangleIndices[3 * i + 2];
}
sphereMesh.startTriangle = 0;
sphereMesh.triangleCount = sphereTriangles.Count();
sphereMesh.startVertex = 0;
sphereMesh.vertexCount = sphereVertices.Count();
vertexCount = sphereVertices.Count();
vertices = sphereVertices.Begin();
triangleCount = sphereTriangles.Count();
triangles = sphereTriangles.Begin();
meshCount = 1;
meshes = &sphereMesh;
sewingLineCount = 0;
sewingLines = nullptr;
}
};
class SoftBody : public Test
{
public:
SoftBody()
{
// Scale and translate the cloth mesh upwards
for (u32 i = 0; i < m_mesh.vertexCount; ++i)
{
m_mesh.vertices[i] *= 2.0f;
m_mesh.vertices[i].y += 10.0f;
}
// Create cloth
b3ClothDef def;
def.mesh = &m_mesh;
def.density = 0.2f;
def.structural = 10000.0f;
def.bending = 0.0f;
def.damping = 0.0f;
m_cloth = new b3Cloth(def);
for (b3Particle* p = m_cloth->GetParticleList().m_head; p; p = p->GetNext())
{
p->SetRadius(0.05f);
p->SetFriction(0.5f);
}
b3Vec3 gravity(0.0f, -9.8f, 0.0f);
m_cloth->SetGravity(gravity);
// Attach a world to the cloth
m_cloth->SetWorld(&m_world);
// Create static shapes
{
b3BodyDef bd;
bd.type = e_staticBody;
b3Body* b = m_world.CreateBody(bd);
static b3BoxHull boxHull(10.0f, 1.0f, 10.0f);
b3HullShape boxShape;
boxShape.m_hull = &boxHull;
b3ShapeDef sd;
sd.shape = &boxShape;
sd.friction = 0.5f;
b->CreateShape(sd);
}
m_clothDragger = new b3ClothDragger(&m_ray, m_cloth);
}
~SoftBody()
{
delete m_clothDragger;
delete m_cloth;
}
// Compute the volume of soft body
float32 ComputeVolume() const
{
const b3ClothMesh* mesh = m_cloth->GetMesh();
b3StackArray<b3Vec3, 256> positions;
positions.Resize(mesh->vertexCount);
for (u32 i = 0; i < mesh->vertexCount; ++i)
{
b3Particle* p = m_cloth->GetVertexParticle(i);
positions[i] = p->GetPosition();
}
b3AABB3 aabb;
aabb.Compute(positions.Begin(), positions.Count());
return aabb.Volume();
}
// Apply pressure forces
// Explanation available in the paper
// "Pressure Model of Soft Body Simulation"
void ApplyPressureForces()
{
const b3ClothMesh* mesh = m_cloth->GetMesh();
// Volume in m^3
float32 V = ComputeVolume();
// Inverse volume
float32 invV = V > 0.0f ? 1.0f / V : 0.0f;
// Apply pressure forces on particles
for (u32 i = 0; i < m_mesh.triangleCount; ++i)
{
u32 i1 = m_mesh.triangles[i].v1;
u32 i2 = m_mesh.triangles[i].v2;
u32 i3 = m_mesh.triangles[i].v3;
b3Particle* p1 = m_cloth->GetVertexParticle(i1);
b3Particle* p2 = m_cloth->GetVertexParticle(i2);
b3Particle* p3 = m_cloth->GetVertexParticle(i3);
b3Vec3 v1 = p1->GetPosition();
b3Vec3 v2 = p2->GetPosition();
b3Vec3 v3 = p3->GetPosition();
b3Vec3 n = b3Cross(v2 - v1, v3 - v1);
// Triangle area
float32 A = n.Normalize();
A *= 0.5f;
// Ideal Gas Approximation
// Number of gas moles
const float32 k_n = 1.0f;
// Ideal Gas Constant [J / K mol]
const float32 k_R = 8.31f;
// Gas temperature in Kelvin
const float32 k_T = 100.0f;
// Pressure in Pascals
float32 P = invV * k_n * k_R * k_T;
// Pressure vector
b3Vec3 Pn = P * n;
// Pressure force
b3Vec3 FP = Pn * A;
// Distribute
p1->ApplyForce(FP);
p2->ApplyForce(FP);
p3->ApplyForce(FP);
}
}
void Step()
{
Test::Step();
ApplyPressureForces();
m_cloth->Step(g_testSettings->inv_hertz);
m_cloth->Draw();
if (m_clothDragger->IsDragging())
{
b3Vec3 pA = m_clothDragger->GetPointA();
b3Vec3 pB = m_clothDragger->GetPointB();
g_draw->DrawPoint(pA, 2.0f, b3Color_green);
g_draw->DrawPoint(pB, 2.0f, b3Color_green);
g_draw->DrawSegment(pA, pB, b3Color_white);
}
extern u32 b3_clothSolverIterations;
g_draw->DrawString(b3Color_white, "Iterations = %d", b3_clothSolverIterations);
float32 E = m_cloth->GetEnergy();
g_draw->DrawString(b3Color_white, "E = %f", E);
}
void MouseMove(const b3Ray3& pw)
{
Test::MouseMove(pw);
if (m_clothDragger->IsDragging() == true)
{
m_clothDragger->Drag();
}
}
void MouseLeftDown(const b3Ray3& pw)
{
Test::MouseLeftDown(pw);
if (m_clothDragger->IsDragging() == false)
{
m_clothDragger->StartDragging();
}
}
void MouseLeftUp(const b3Ray3& pw)
{
Test::MouseLeftUp(pw);
if (m_clothDragger->IsDragging() == true)
{
m_clothDragger->StopDragging();
}
}
static Test* Create()
{
return new SoftBody();
}
b3SphereClothMesh m_mesh;
b3Cloth* m_cloth;
b3ClothDragger* m_clothDragger;
};
#endif