Skip to content

Commit

Permalink
updated version
Browse files Browse the repository at this point in the history
  • Loading branch information
Eftychios Sifakis committed May 5, 2022
1 parent 331497a commit 82dc905
Show file tree
Hide file tree
Showing 3 changed files with 39 additions and 81 deletions.
4 changes: 2 additions & 2 deletions Collision-Tests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,11 @@ set (CMAKE_CXX_STANDARD 14)

project (Collision-Tests)

set (CMAKE_CXX_FLAGS "-Wno-deprecated ${CMAKE_CXX_FLAGS}")
set (CMAKE_CXX_FLAGS "-fopenmp -Wno-deprecated ${CMAKE_CXX_FLAGS}")

add_definitions(-DBUILD_COMPONENT_SRC_PREFIX="" -DBUILD_OPTLEVEL_DEV)

set(USD_ROOT_DIR ~/USD)
set(USD_ROOT_DIR /nobackup/sifakis/opt)
set(USD_LIB_DIR ${USD_ROOT_DIR}/lib/)
set(USD_INC_DIR ${USD_ROOT_DIR}/include/)

Expand Down
114 changes: 36 additions & 78 deletions Collision-Tests/tests/GroundCollision3D/FiniteElementMesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@
#include <Eigen/Dense>

#define USE_LINEAR_ELASTICITY
//#define USE_ST_VENANT_KIRCHHOFF
//#define USE_COROTATED_ELASTICITY

template<class T>
Expand Down Expand Up @@ -37,7 +36,10 @@ struct FiniteElementMesh : public AnimatedTetrahedronMesh<T>
std::vector<Vector3> m_particleV;
std::vector<Matrix33> m_DmInverse;
std::vector<T> m_restVolume;


mutable std::vector<Matrix33> Ue,Ve;
mutable std::vector<Vector3> vSigmae;

std::vector<int> m_surfaceParticles;
mutable std::vector<bool> m_collisionActive;
mutable std::vector<Vector3> m_collisionTarget;
Expand Down Expand Up @@ -70,6 +72,11 @@ struct FiniteElementMesh : public AnimatedTetrahedronMesh<T>

void addElasticForce(std::vector<Vector3>& f) const
{
Ue.resize(m_meshElements.size());
Ve.resize(m_meshElements.size());
vSigmae.resize(m_meshElements.size());

#pragma omp parallel for
for(int e = 0; e < m_meshElements.size(); e++)
{
const auto& element = m_meshElements[e];
Expand All @@ -82,9 +89,12 @@ struct FiniteElementMesh : public AnimatedTetrahedronMesh<T>

// Compute SVD
Eigen::JacobiSVD<Matrix33> svd(F, Eigen::ComputeFullU | Eigen::ComputeFullV);
Matrix33 U = svd.matrixU();
Matrix33 V = svd.matrixV();
Vector3 vSigma = svd.singularValues();
Matrix33& U = Ue[e];
Matrix33& V = Ve[e];
Vector3& vSigma = vSigmae[e];
U = svd.matrixU();
V = svd.matrixV();
vSigma = svd.singularValues();
if ( U.determinant() < 0. ) {
if ( V.determinant() < 0. ) {
// Both determinants negative, just negate last column on both
Expand All @@ -108,22 +118,28 @@ struct FiniteElementMesh : public AnimatedTetrahedronMesh<T>
for (int v = 0; v < 3; v++)
vSigma[v] = std::max<T>(m_singularValueThreshold, vSigma[v]);
Matrix33 Sigma = vSigma.asDiagonal();
F = U * Sigma * V.transpose();
}

for(int e = 0; e < m_meshElements.size(); e++)
{
const auto& element = m_meshElements[e];

// Use precomputed SVD
const Matrix33& U = Ue[e];
const Matrix33& V = Ue[e];
const Vector3& vSigma = vSigmae[e];
const Matrix33 F = U * vSigma.asDiagonal() * V.transpose();

#ifdef USE_LINEAR_ELASTICITY
Matrix33 strain = .5 * (F + F.transpose()) - Matrix33::Identity();
Matrix33 P = 2. * m_mu * strain + m_lambda * strain.trace() * Matrix33::Identity();
#endif

#ifdef USE_ST_VENANT_KIRCHHOFF
Matrix22 E = .5 * ( F.transpose() * F - Matrix22::Identity());
Matrix22 P = F * (2. * m_mu * E + m_lambda * E.trace() * Matrix22::Identity());
#endif

#ifdef USE_COROTATED_ELASTICITY
Vector2 vStrain = vSigma - Vector2::Ones();
Vector2 vP = 2. * m_mu * vStrain + m_lambda * vStrain.sum() * Vector2::Ones();
Matrix22 P = U * vP.asDiagonal() * V.transpose();
Vector3 vStrain = vSigma - Vector3::Ones();
Vector3 vP = 2. * m_mu * vStrain + m_lambda * vStrain.sum() * Vector3::Ones();
Matrix33 P = U * vP.asDiagonal() * V.transpose();
#endif

Matrix33 H = -m_restVolume[e] * P * m_DmInverse[e].transpose();
Expand Down Expand Up @@ -175,41 +191,10 @@ struct FiniteElementMesh : public AnimatedTetrahedronMesh<T>
{
const auto& element = m_meshElements[e];

// Compute deformation gradient
Matrix33 Ds;
for(int j = 0; j < 3; j++)
Ds.col(j) = m_particleX[element[j+1]]-m_particleX[element[0]];
Matrix33 F = Ds * m_DmInverse[e];

// Compute SVD
Eigen::JacobiSVD<Matrix33> svd(F, Eigen::ComputeFullU | Eigen::ComputeFullV);
Matrix33 U = svd.matrixU();
Matrix33 V = svd.matrixV();
Vector3 vSigma = svd.singularValues();
if ( U.determinant() < 0. ) {
if ( V.determinant() < 0. ) {
// Both determinants negative, just negate last column on both
U.col(2) *= -1.f;
V.col(2) *= -1.f;
}
else {
// Only U has negative determinant, negate last column and second singular value
U.col(2) *= -1.f;
vSigma[2] = -vSigma[2];
}
}
else
if ( V.determinant() < 0.) {
// Only V has negative determinant, negate last column and second singular value
V.col(2) *= -1.f;
vSigma[2] = -vSigma[2];
}

// Apply thresholding of singular values, and re-constitute F
for (int v = 0; v < 3; v++)
vSigma[v] = std::max<T>(m_singularValueThreshold, vSigma[v]);
Matrix33 Sigma = vSigma.asDiagonal();
F = U * Sigma * V.transpose();
// Use precomputed SVD
const Matrix33& U = Ue[e];
const Matrix33& V = Ue[e];
const Vector3& vSigma = vSigmae[e];

// Compute differential(s)
Matrix33 dDs;
Expand All @@ -222,38 +207,11 @@ struct FiniteElementMesh : public AnimatedTetrahedronMesh<T>
Matrix33 dP = 2. * m_mu * dstrain + m_lambda * dstrain.trace() * Matrix33::Identity();
#endif

#ifdef USE_ST_VENANT_KIRCHHOFF
Matrix22 E = .5 * ( F.transpose() * F - Matrix22::Identity());
Matrix22 dE = .5 * ( dF.transpose() * F + F.transpose() * dF);
Matrix22 dP = dF * (2. * m_mu * E + m_lambda * E.trace() * Matrix22::Identity()) +
F * (2. * m_mu * dE + m_lambda * dE.trace() * Matrix22::Identity());

#endif

#ifdef USE_COROTATED_ELASTICITY
// Construct diagonalized dP/dF tensor
DiagonalizedStressDerivative dsd;
dsd.A = Matrix22::Ones() * m_lambda + Matrix22::Identity() * (2. * m_mu);

Vector2 vStrain = vSigma - Vector2::Ones();
T q = ( m_mu + m_lambda) * std::max<T>( vStrain.sum(), 0. ) / ( vSigma.sum() + 1e-4 ); // Definiteness fix
dsd.B12(0, 0) = dsd.B12(1, 1) = m_mu + q;
dsd.B12(0, 1) = dsd.B12(1, 0) = m_mu - q;

// Apply tensor (with required rotations)
Matrix22 dF_hat = U.transpose() * dF * V;

Vector2 vdP_A = dsd.A * Vector2( dF_hat(0, 0), dF_hat(1, 1) );
Vector2 vdP_B12 = dsd.B12 * Vector2( dF_hat(0, 1), dF_hat(1, 0) );
Matrix22 dP_hat;
dP_hat(0, 0) = vdP_A[0];
dP_hat(1, 1) = vdP_A[1];
dP_hat(0, 1) = vdP_B12[0];
dP_hat(1, 0) = vdP_B12[1];

Matrix22 dP = U * dP_hat * V.transpose();
// Implementing stiffness via projective dynamics instead!
Matrix33 dP = 2. * m_mu * dF;
#endif

Matrix33 dH = (scale * m_restVolume[e]) * dP * m_DmInverse[e].transpose();

for(int j = 0; j < 3; j++){
Expand Down
2 changes: 1 addition & 1 deletion Collision-Tests/tests/GroundCollision3D/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@ int main(int argc, char *argv[])
LatticeMesh<float> simulationMesh;
simulationMesh.m_cellSize = { 10, 10, 10 };
simulationMesh.m_gridDX = 0.1;
simulationMesh.m_nFrames = 50;
simulationMesh.m_nFrames = 10;
simulationMesh.m_subSteps = 1;
simulationMesh.m_frameDt = 0.04;

Expand Down

0 comments on commit 82dc905

Please sign in to comment.