A new PBD FEMTriangleConstraint

Interactive Medical Simulation Toolkit (imstk)
1

I want to integrate a new constraint.But I do not know how to compute in XPBD mode.
in imstk the FEMTetConstraint compute:

m_compliance = 1.0 / (config->m_lambda + 2 * config->m_mu);

but in PBD FEMTriangleConstraint, is a Orthotropic material by these parameter:

const Real youngsModulusX
const Real youngsModulusY
const Real youngsModulusShear
const Real poissonRatioXY
const Real poissonRatioYX
how to compute m_compliance .

I get the Constraint form here:FEMTriangleConstraint

this is the class and slove:

class FEMTriangleConstraint : public Constraint
{
public:
static int TYPE_ID;
Real m_area;
Matrix2r m_invRestMat;

  FEMTriangleConstraint() : Constraint(3) {}
  virtual int &getTypeId() const { return TYPE_ID; }

  virtual bool initConstraint(SimulationModel &model, const unsigned int particle1, const unsigned int particle2,
  	const unsigned int particle3);
  virtual bool solvePositionConstraint(SimulationModel &model, const unsigned int iter);

};

bool FEMTriangleConstraint::initConstraint(SimulationModel &model, const unsigned int particle1, const unsigned int particle2,
const unsigned int particle3)
{
m_bodies[0] = particle1;
m_bodies[1] = particle2;
m_bodies[2] = particle3;

ParticleData &pd = model.getParticles();

Vector3r &x1 = pd.getPosition0(particle1);
Vector3r &x2 = pd.getPosition0(particle2);
Vector3r &x3 = pd.getPosition0(particle3);

return PositionBasedDynamics::init_FEMTriangleConstraint(x1, x2, x3, m_area, m_invRestMat);
}

bool FEMTriangleConstraint::solvePositionConstraint(SimulationModel &model, const unsigned int iter)
{
ParticleData &pd = model.getParticles();

const unsigned i1 = m_bodies[0];
const unsigned i2 = m_bodies[1];
const unsigned i3 = m_bodies[2];

Vector3r &x1 = pd.getPosition(i1);
Vector3r &x2 = pd.getPosition(i2);
Vector3r &x3 = pd.getPosition(i3);

const Real invMass1 = pd.getInvMass(i1);
const Real invMass2 = pd.getInvMass(i2);
const Real invMass3 = pd.getInvMass(i3);

Vector3r corr1, corr2, corr3;
const bool res = PositionBasedDynamics::solve_FEMTriangleConstraint(
x1, invMass1,
x2, invMass2,
x3, invMass3,
m_area,
m_invRestMat,
model.getValue(SimulationModel::CLOTH_STIFFNESS_XX),
model.getValue(SimulationModel::CLOTH_STIFFNESS_YY),
model.getValue(SimulationModel::CLOTH_STIFFNESS_XY),
model.getValue(SimulationModel::CLOTH_POISSON_RATIO_XY),
model.getValue(SimulationModel::CLOTH_POISSON_RATIO_YX),
corr1, corr2, corr3);

if (res)
{
if (invMass1 != 0.0)
x1 += corr1;
if (invMass2 != 0.0)
x2 += corr2;
if (invMass3 != 0.0)
x3 += corr3;
}
return res;
}

// ----------------------------------------------------------------------------------------------
bool PositionBasedDynamics::init_FEMTriangleConstraint(
const Vector3r &p0,
const Vector3r &p1,
const Vector3r &p2,
Real &area,
Matrix2r &invRestMat)
{
Vector3r normal0 = (p1 - p0).cross(p2 - p0);
area = normal0.norm() * static_cast(0.5);

Vector3r axis0_1 = p1 - p0;
axis0_1.normalize();
Vector3r axis0_2 = normal0.cross(axis0_1);
axis0_2.normalize();

Vector2r p[3];
p[0] = Vector2r(p0.dot(axis0_2), p0.dot(axis0_1));
p[1] = Vector2r(p1.dot(axis0_2), p1.dot(axis0_1));
p[2] = Vector2r(p2.dot(axis0_2), p2.dot(axis0_1));

Matrix2r P;
P(0, 0) = p[0][0] - p[2][0];
P(1, 0) = p[0][1] - p[2][1];
P(0, 1) = p[1][0] - p[2][0];
P(1, 1) = p[1][1] - p[2][1];

const Real det = P.determinant();
if (fabs(det) > eps)
{
invRestMat = P.inverse();
return true;
}
return false;
}

// ----------------------------------------------------------------------------------------------
bool PositionBasedDynamics::solve_FEMTriangleConstraint(
const Vector3r &p0, Real invMass0,
const Vector3r &p1, Real invMass1,
const Vector3r &p2, Real invMass2,
const Real &area,
const Matrix2r &invRestMat,
const Real youngsModulusX,
const Real youngsModulusY,
const Real youngsModulusShear,
const Real poissonRatioXY,
const Real poissonRatioYX,
Vector3r &corr0, Vector3r &corr1, Vector3r &corr2)
{
// Orthotropic elasticity tensor
Matrix3r C;
C.setZero();
C(0, 0) = youngsModulusX / (static_cast(1.0) - poissonRatioXYpoissonRatioYX);
C(0, 1) = youngsModulusXpoissonRatioYX / (static_cast(1.0) - poissonRatioXYpoissonRatioYX);
C(1, 1) = youngsModulusY / (static_cast(1.0) - poissonRatioXYpoissonRatioYX);
C(1, 0) = youngsModulusYpoissonRatioXY / (static_cast(1.0) - poissonRatioXYpoissonRatioYX);
C(2, 2) = youngsModulusShear;

// Determine \partial x/\partial m_i
Eigen::Matrix<Real, 3, 2> F;
const Vector3r p13 = p0 - p2;
const Vector3r p23 = p1 - p2;
F(0,0) = p13[0] * invRestMat(0,0) + p23[0] * invRestMat(1,0);
F(0,1) = p13[0] * invRestMat(0,1) + p23[0] * invRestMat(1,1);
F(1,0) = p13[1] * invRestMat(0,0) + p23[1] * invRestMat(1,0);
F(1,1) = p13[1] * invRestMat(0,1) + p23[1] * invRestMat(1,1);
F(2,0) = p13[2] * invRestMat(0,0) + p23[2] * invRestMat(1,0);
F(2,1) = p13[2] * invRestMat(0,1) + p23[2] * invRestMat(1,1);

// epsilon = 0.5(F^T * F - I)
Matrix2r epsilon;
epsilon(0,0) = static_cast(0.5)(F(0,0) * F(0,0) + F(1,0) * F(1,0) + F(2,0) * F(2,0) - static_cast(1.0)); // xx
epsilon(1,1) = static_cast(0.5)(F(0,1) * F(0,1) + F(1,1) * F(1,1) + F(2,1) * F(2,1) - static_cast(1.0)); // yy
epsilon(0,1) = static_cast(0.5)*(F(0,0) * F(0,1) + F(1,0) * F(1,1) + F(2,0) * F(2,1)); // xy
epsilon(1,0) = epsilon(0,1);

// P(F) = det(F) * C*E * F^-T => E = green strain
Matrix2r stress;
stress(0,0) = C(0,0) * epsilon(0,0) + C(0,1) * epsilon(1,1) + C(0,2) * epsilon(0,1);
stress(1,1) = C(1,0) * epsilon(0,0) + C(1,1) * epsilon(1,1) + C(1,2) * epsilon(0,1);
stress(0,1) = C(2,0) * epsilon(0,0) + C(2,1) * epsilon(1,1) + C(2,2) * epsilon(0,1);
stress(1,0) = stress(0,1);

const Eigen::Matrix<Real, 3, 2> piolaKirchhoffStres = F * stress;

Real psi = 0.0;
for (unsigned char j = 0; j < 2; j++)
for (unsigned char k = 0; k < 2; k++)
psi += epsilon(j,k) * stress(j,k);
psi = static_cast(0.5)psi;
Real energy = areapsi;

// compute gradient
Eigen::Matrix<Real, 3, 2> H = area * piolaKirchhoffStres * invRestMat.transpose();
Vector3r gradC[3];
for (unsigned char j = 0; j < 3; ++j)
{
gradC[0][j] = H(j,0);
gradC[1][j] = H(j,1);
}
gradC[2] = -gradC[0] - gradC[1];

Real sum_normGradC = invMass0 * gradC[0].squaredNorm();
sum_normGradC += invMass1 * gradC[1].squaredNorm();
sum_normGradC += invMass2 * gradC[2].squaredNorm();

// exit early if required
if (fabs(sum_normGradC) > eps)
{
// compute scaling factor
const Real s = energy / sum_normGradC;

  // update positions
  corr0 = -(s*invMass0) * gradC[0];
  corr1 = -(s*invMass1) * gradC[1];
  corr2 = -(s*invMass2) * gradC[2];

  return true;

}

return false;
}

1 Like
2

We haven’t specifically looked into the compliance matrix of the orthotropic material specifically but based on the literature it is the inverse of the stiffness matrix (see xPBD paper). Also, a quick google search gave us this.

In imstk, PBD and xPBD share the same framework and the implementation and the only thing that changes is the solver. To implement a new xpbd constraint, one has to

  1. implement a virtual function PbdConstraint::computeValueAndGradient.
  2. specify a value for m_compliance

Currently, we use a scalar for m_compliance for simplicity but we are planning to specify constraint specific compliance matrix similar to Eq (28) in the xpbd paper.

Hope this helps. Please let us know if we can be of further assistance in implementation or in understanding the literature. We will also be happy to review your code.

Best.