Revision control
++ // Copyright (c) 2001-2004 NovodeX AG. All rights reserved.
++ // Copyright (c) 2004-2008 AGEIA Technologies, Inc. All rights reserved.
++ // SPDX-FileCopyrightText: Copyright (c) 2008-2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
++ // SPDX-License-Identifier: Apache-2.0
++
++ #ifndef GU_SLOW_GJK_EPA_H
++ #define GU_SLOW_GJK_EPA_H
++
++ #include "foundation/PxBasicTemplates.h"
++ #include "foundation/PxVec3.h"
++ #include "foundation/PxVec4.h"
++ #include "GuBounds.h"
++
++ namespace physx
++ {
++ namespace Gu
++ {
++ // The reference implementation of GJK-EPA algorithm
++ namespace RefGjkEpa
++ {
++ template <typename Support>
++ struct Convex
++ {
++ static const PxU32 MAX_VERTS = 32;
++
++ PX_CUDA_CALLABLE
++ Convex(const Support& s, const PxTransform& p)
++ :
++ mS(s), mPose(p), mNumVerts(0)
++ {
++ const PxVec3 X(1, 0, 0), Y(0, 1, 0), Z(0, 0, 1);
++ const PxBounds3 bounds(PxVec3(mS.supportLocal(-X).x, mS.supportLocal(-Y).y, mS.supportLocal(-Z).z),
++ PxVec3(mS.supportLocal(X).x, mS.supportLocal(Y).y, mS.supportLocal(Z).z));
++ PxReal maxExtent = bounds.getDimensions().maxElement();
++ mAccuracy = PxMax(maxExtent * 0.01f, FLT_EPSILON);
++ }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxReal getAccuracy() const
++ { return mAccuracy; }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxBounds3 getBounds() const
++ {
++ const PxVec3 X(mPose.q.getInvBasisVector0()),
++ Y(mPose.q.getInvBasisVector1()),
++ Z(mPose.q.getInvBasisVector2());
++
++ return PxBounds3(PxVec3(mPose.transform(mS.supportLocal(-X)).x,
++ mPose.transform(mS.supportLocal(-Y)).y,
++ mPose.transform(mS.supportLocal(-Z)).z),
++ PxVec3(mPose.transform(mS.supportLocal(X)).x,
++ mPose.transform(mS.supportLocal(Y)).y,
++ mPose.transform(mS.supportLocal(Z)).z));
++ }
++
++ PX_CUDA_CALLABLE
++ PxU8 supportIndex(const PxVec3& dir)
++ {
++ PxVec3 d = dir;
++ PxVec3 v = mPose.transform(mS.supportLocal(mPose.rotateInv(d)));
++ PxU32 closest = MAX_VERTS;
++ const PxReal WELD_DIST_SQ = mAccuracy * mAccuracy;
++ PxReal distSq = FLT_MAX;
++ for (PxI32 i = mNumVerts - 1; i >= 0; --i)
++ {
++ PxReal dSq = (mVerts[i] - v).magnitudeSquared();
++
++ if (dSq < WELD_DIST_SQ)
++ return PxU8(i);
++
++ if (distSq > dSq)
++ {
++ distSq = dSq;
++ closest = i;
++ }
++ }
++
++ PX_ASSERT(closest < MAX_VERTS || mNumVerts == 0);
++
++ if (mNumVerts == MAX_VERTS)
++ return PxU8(closest);
++
++ mVerts[mNumVerts++] = v;
++
++ return PxU8(mNumVerts - 1);
++ }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ const PxVec3& supportVertex(PxU8 index) const
++ { PX_ASSERT(index < mNumVerts); return mVerts[index]; }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ const PxTransform& getPose() const
++ { return mPose; }
++
++ private:
++
++ const Support& mS;
++ const PxTransform& mPose;
++ PxVec3 mVerts[MAX_VERTS];
++ PxU32 mNumVerts;
++ PxF32 mAccuracy;
++ };
++
++ template <typename Support>
++ struct GjkDistance
++ {
++ PX_CUDA_CALLABLE
++ GjkDistance(Convex<Support>& convexA, Convex<Support>& convexB)
++ :
++ mConvexA(convexA), mConvexB(convexB),
++ mNumVerts(0), mNumBestVerts(0), mIteration(0), mClosest(FLT_MAX)
++ {
++ mVerts[0] = mVerts[1] = mVerts[2] = mVerts[3] = Vert::make(PxVec3(0), 0, 0, 0);
++ PxVec3 dir = mConvexB.getBounds().getCenter() - mConvexA.getBounds().getCenter();
++ if (dir.normalize() < FLT_EPSILON) dir = PxVec3(1, 0, 0);
++ addPoint(dir);
++ }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxReal getAccuracy() const
++ { return PxMin(mConvexA.getAccuracy(), mConvexB.getAccuracy()); }
++
++ PX_CUDA_CALLABLE
++ // adds a new point (the difference of the shapes support points in a given direction)
++ // to the simplex. returns the projection of the difference of the support points on
++ // the direction. positive values mean a gap, negative - an overlap. FLT_MAX signals
++ // to finalize the algorithm. mIteration > MAX_ITERATIONS indicate that the algorithm
++ // got stuck. most likely adding several same points over and over, but unable to
++ // finish because of the floating point precision problems
++ PxReal addPoint(const PxVec3& dir)
++ {
++ if (++mIteration > MAX_ITERATIONS)
++ return FLT_MAX;
++
++ const PxU8 aI = mConvexA.supportIndex(dir), bI = mConvexB.supportIndex(-dir);
++ const PxVec3 p = mConvexA.supportVertex(aI) - mConvexB.supportVertex(bI);
++
++ for (PxU32 i = 0; i < mNumVerts; ++i)
++ {
++ const Vert& v = mVerts[i];
++ if (v.aI == aI && v.bI == bI)
++ return FLT_MAX;
++ }
++
++ mVerts[mNumVerts++] = Vert::make(p, 0.0f, aI, bI);
++
++ return p.dot(-dir);
++ }
++
++ PX_CUDA_CALLABLE
++ void computeClosest()
++ {
++ Vert verts[4]; verts[0] = mVerts[0]; verts[1] = mVerts[1]; verts[2] = mVerts[2]; verts[3] = mVerts[3];
++ verts[0].s = verts[1].s = verts[2].s = verts[3].s = 0;
++ PxU32 numVerts = mNumVerts;
++
++ switch (numVerts)
++ {
++ case 1:
++ {
++ verts[0].s = 1.0f;
++ break;
++ }
++ case 2:
++ {
++ const PxVec3 a = verts[0].p, b = verts[1].p;
++ const PxVec3 ba = a - b;
++ const PxReal sba = ba.dot(-b);
++ if (sba <= 0)
++ {
++ verts[0] = verts[1];
++ verts[0].s = 1.0f;
++ numVerts = 1;
++ break;
++ }
++ verts[0].s = sba / ba.magnitudeSquared();
++ verts[1].s = 1.0f - verts[0].s;
++ break;
++ }
++ case 3:
++ {
++ const PxVec3 a = verts[0].p, b = verts[1].p, c = verts[2].p;
++ const PxVec3 ca = a - c, cb = b - c;
++ const PxReal sca = ca.dot(-c), scb = cb.dot(-c);
++ if (sca <= 0 && scb <= 0)
++ {
++ verts[0] = verts[2];
++ verts[0].s = 1.0f;
++ numVerts = 1;
++ break;
++ }
++ const PxVec3 abc = (b - a).cross(c - a);
++ const PxReal iabc = 1.0f / abc.magnitudeSquared();
++ const PxVec3 pabc = abc * abc.dot(a) * iabc;
++ const PxReal tca = abc.dot((c - pabc).cross(a - pabc));
++ if (sca > 0 && tca <= 0)
++ {
++ verts[1] = verts[2];
++ verts[0].s = sca / ca.magnitudeSquared();
++ verts[1].s = 1.0f - verts[0].s;
++ numVerts = 2;
++ break;
++ }
++ const PxReal tbc = abc.dot((b - pabc).cross(c - pabc));
++ if (scb > 0 && tbc <= 0)
++ {
++ verts[0] = verts[1];
++ verts[1] = verts[2];
++ verts[0].s = scb / cb.magnitudeSquared();
++ verts[1].s = 1.0f - verts[0].s;
++ numVerts = 2;
++ break;
++ }
++ verts[0].s = tbc * iabc;
++ verts[1].s = tca * iabc;
++ verts[2].s = 1.0f - verts[0].s - verts[1].s;
++ if (abc.dot(a) > 0) PxSwap(verts[0], verts[1]);
++ break;
++ }
++ case 4:
++ {
++ const PxVec3 a = verts[0].p, b = verts[1].p, c = verts[2].p, d = verts[3].p;
++ const PxVec3 da = a - d, db = b - d, dc = c - d;
++ const PxReal sda = da.dot(-d), sdb = db.dot(-d), sdc = dc.dot(-d);
++ if (sda <= 0 && sdb <= 0 && sdc <= 0)
++ {
++ verts[0] = verts[3];
++ verts[0].s = 1.0f;
++ numVerts = 1;
++ break;
++ }
++ const PxVec3 dab = (a - d).cross(b - d);
++ const PxReal idab = 1.0f / dab.magnitudeSquared();
++ const PxVec3 pdab = dab * dab.dot(d) * idab;
++ const PxVec3 dbc = (b - d).cross(c - d);
++ const PxReal idbc = 1.0f / dbc.magnitudeSquared();
++ const PxVec3 pdbc = dbc * dbc.dot(d) * idbc;
++ const PxVec3 dca = (c - d).cross(a - d);
++ const PxReal idca = 1.0f / dca.magnitudeSquared();
++ const PxVec3 pdca = dca * dca.dot(d) * idca;
++ const PxReal tda = dab.dot((d - pdab).cross(a - pdab));
++ const PxReal tad = dca.dot((a - pdca).cross(d - pdca));
++ if (sda > 0 && tda <= 0 && tad <= 0)
++ {
++ verts[1] = verts[3];
++ verts[0].s = sda / da.magnitudeSquared();
++ verts[1].s = 1.0f - verts[0].s;
++ numVerts = 2;
++ break;
++ }
++ const PxReal tdb = dbc.dot((d - pdbc).cross(b - pdbc));
++ const PxReal tbd = dab.dot((b - pdab).cross(d - pdab));
++ if (sdb > 0 && tdb <= 0 && tbd <= 0)
++ {
++ verts[0] = verts[1];
++ verts[1] = verts[3];
++ verts[0].s = sdb / db.magnitudeSquared();
++ verts[1].s = 1.0f - verts[0].s;
++ numVerts = 2;
++ break;
++ }
++ const PxReal tcd = dbc.dot((c - pdbc).cross(d - pdbc));
++ const PxReal tdc = dca.dot((d - pdca).cross(c - pdca));
++ if (sdc > 0 && tdc <= 0 && tcd <= 0)
++ {
++ verts[0] = verts[2];
++ verts[1] = verts[3];
++ verts[0].s = sdc / dc.magnitudeSquared();
++ verts[1].s = 1.0f - verts[0].s;
++ numVerts = 2;
++ break;
++ }
++ const PxVec3 abc = (b - a).cross(c - a);
++ if (tda > 0 && tbd > 0 && dab.dot(d) < 0 && abc.dot(dab) > 0)
++ {
++ verts[2] = verts[3];
++ verts[0].s = tbd * idab;
++ verts[1].s = tda * idab;
++ verts[2].s = 1.0f - verts[0].s - verts[1].s;
++ numVerts = 3;
++ break;
++ }
++ if (tdb > 0 && tcd > 0 && dbc.dot(d) < 0 && abc.dot(dbc) > 0)
++ {
++ verts[0] = verts[1];
++ verts[1] = verts[2];
++ verts[2] = verts[3];
++ verts[0].s = tcd * idbc;
++ verts[1].s = tdb * idbc;
++ verts[2].s = 1.0f - verts[0].s - verts[1].s;
++ numVerts = 3;
++ break;
++ }
++ if (tdc > 0 && tad > 0 && dca.dot(d) < 0 && abc.dot(dca) > 0)
++ {
++ verts[1] = verts[3];
++ verts[0].s = tdc * idca;
++ verts[2].s = tad * idca;
++ verts[1].s = 1.0f - verts[0].s - verts[2].s;
++ numVerts = 3;
++ break;
++ }
++ break;
++ }
++ }
++
++ PxVec3 closest = verts[0].p * verts[0].s + verts[1].p * verts[1].s + verts[2].p * verts[2].s + verts[3].p * verts[3].s;
++
++ if (closest.magnitude() < mClosest.magnitude() - FLT_EPSILON)
++ {
++ mClosest = closest;
++ mVerts[0] = verts[0]; mVerts[1] = verts[1]; mVerts[2] = verts[2]; mVerts[3] = verts[3];
++ mNumVerts = numVerts;
++ }
++ else
++ --mNumVerts;
++ }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxReal getDist() const
++ { return mClosest.magnitude(); }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxVec3 getDir() const
++ { return -mClosest.getNormalized(); }
++
++ PX_CUDA_CALLABLE
++ void computePoints(PxVec3& pointA, PxVec3& pointB)
++ {
++ PxVec3 pA(0), pB(0);
++ for (PxU32 i = 0; i < mNumVerts; ++i)
++ {
++ pA += mConvexA.supportVertex(mVerts[i].aI) * mVerts[i].s;
++ pB += mConvexB.supportVertex(mVerts[i].bI) * mVerts[i].s;
++ }
++ pointA = pA; pointB = pB;
++ }
++
++ private:
++
++ static const PxU32 MAX_ITERATIONS = 100;
++
++ struct Vert
++ {
++ PxVec3 p;
++ PxReal s;
++ PxU8 aI, bI;
++
++ PX_CUDA_CALLABLE PX_INLINE
++ static Vert make(const PxVec3& p, PxReal s, PxU8 aI, PxU8 bI)
++ { Vert v; v.p = p; v.s = s; v.aI = aI; v.bI = bI; return v; }
++ };
++
++ Convex<Support>& mConvexA;
++ Convex<Support>& mConvexB;
++ PxU32 mNumVerts, mNumBestVerts, mIteration;
++ Vert mVerts[4];
++ PxVec3 mClosest;
++ };
++
++ template <typename Support>
++ PX_CUDA_CALLABLE
++ static PxReal computeGjkDistance(const Support& a, const Support& b, const PxTransform& poseA, const PxTransform& poseB,
++ PxReal maxDist, PxVec3& pointA, PxVec3& pointB, PxVec3& axis)
++ {
++ Convex<Support> convexA(a, poseA), convexB(b, poseB);
++ GjkDistance<Support> gjk(convexA, convexB);
++
++ const PxReal distEps = 0.01f * gjk.getAccuracy();
++
++ while (true)
++ {
++ gjk.computeClosest();
++
++ const PxReal closestDist = gjk.getDist();
++ if (closestDist < distEps)
++ return 0;
++
++ const PxVec3 dir = gjk.getDir();
++ const PxReal proj = gjk.addPoint(dir);
++
++ if (proj >= closestDist - distEps)
++ {
++ PxVec3 pA, pB;
++ gjk.computePoints(pA, pB);
++ pointA = pA; pointB = pB;
++ axis = -dir;
++ return closestDist;
++ }
++
++ if (proj > maxDist)
++ return FLT_MAX;
++ }
++ }
++
++ template <typename Support>
++ struct EpaDepth
++ {
++ PX_CUDA_CALLABLE
++ EpaDepth(Convex<Support>& convexA, Convex<Support>& convexB)
++ :
++ mConvexA(convexA), mConvexB(convexB),
++ mNumVerts(0), mNumFaces(0),
++ mClosest(Face::make(PxVec4(PxZero), 0, 0, 0)),
++ mProj(-FLT_MAX)
++ {
++ PxVec3 dir0 = mConvexB.getBounds().getCenter() - mConvexA.getBounds().getCenter();
++ if (dir0.normalize() < FLT_EPSILON) dir0 = PxVec3(1, 0, 0);
++ PxVec3 dir1, dir2; PxComputeBasisVectors(dir0, dir1, dir2);
++ addPoint(dir0); addPoint(-dir0); addPoint(dir1);
++ if (mNumVerts < 3) addPoint(-dir1);
++ if (mNumVerts < 3) addPoint(dir2);
++ if (mNumVerts < 3) addPoint(-dir2);
++ mProj = -FLT_MAX;
++ }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ bool isValid() const
++ { return mNumFaces > 0; }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxReal getAccuracy() const
++ { return PxMin(mConvexA.getAccuracy(), mConvexB.getAccuracy()); }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxVec3 supportVertex(PxU8 aI, PxU8 bI)
++ { return mConvexA.supportVertex(aI) - mConvexB.supportVertex(bI); }
++
++ PX_CUDA_CALLABLE
++ bool makePlane(PxU8 i0, PxU8 i1, PxU8 i2, PxVec4& plane)
++ {
++ const Vert v0 = mVerts[i0], v1 = mVerts[i1], v2 = mVerts[i2];
++ const PxVec3 a = supportVertex(v0.aI, v0.bI), b = supportVertex(v1.aI, v1.bI), c = supportVertex(v2.aI, v2.bI);
++ const PxVec3 abc = (b - a).cross(c - a) + (c - b).cross(a - b) + (a - c).cross(b - c);
++ const PxReal labc = abc.magnitude();
++
++ const PxReal normalEps = FLT_EPSILON;
++ if (labc < normalEps)
++ return false;
++
++ const PxVec3 n = abc / labc;
++ const PxReal d = n.dot(-a);
++
++ plane = PxVec4(n, d);
++
++ return true;
++ }
++
++ PX_CUDA_CALLABLE
++ // adds a new point (the difference of the shapes support points in a given direction)
++ // to the polytope. removes all faces below the new point, keeps the track of open edges
++ // created by the face removal, and then creates a fan of new faces each containing the
++ // new point and one of the open edges. returns the projection of the difference of the
++ // support points on the direction. FLT_MAX signals to finalize the algorithm.
++ PxReal addPoint(const PxVec3& dir)
++ {
++ if (mNumVerts > MAX_VERTS - 1)
++ return FLT_MAX;
++
++ const PxU8 aI = mConvexA.supportIndex(dir), bI = mConvexB.supportIndex(-dir);
++ const PxVec3 p = mConvexA.supportVertex(aI) - mConvexB.supportVertex(bI);
++
++ PxReal proj = p.dot(-dir);
++ if (proj > mProj)
++ {
++ mProj = proj;
++ mBest = mClosest;
++ }
++
++ for (PxU32 i = 0; i < mNumVerts; ++i)
++ if (aI == mVerts[i].aI && bI == mVerts[i].bI)
++ return FLT_MAX;
++
++ Vert& v = mVerts[mNumVerts++];
++ v = Vert::make(aI, bI);
++
++ if (mNumVerts < 3)
++ return 0;
++
++ if (mNumVerts == 3)
++ {
++ PxVec4 plane0, plane1;
++ if (!makePlane(0, 1, 2, plane0) || !makePlane(1, 0, 2, plane1))
++ return FLT_MAX;
++
++ mFaces[mNumFaces++] = Face::make(plane0, 0, 1, 2);
++ mFaces[mNumFaces++] = Face::make(plane1, 1, 0, 2);
++
++ return 0;
++ }
++
++ const PxReal testEps = 0.01f * getAccuracy();
++ Edge edges[MAX_EDGES];
++ PxU32 numEdges = 0;
++ for (PxI32 i = mNumFaces - 1; i >= 0; --i)
++ {
++ Face& f = mFaces[i];
++ if (f.plane.dot(PxVec4(p, 1)) > testEps)
++ {
++ for (PxU32 j = 0; j < 3; ++j)
++ {
++ bool add = true;
++ const PxU8 v0 = f.v[j], v1 = f.v[(j + 1) % 3];
++ for (PxI32 k = numEdges - 1; k >= 0; --k)
++ {
++ const Edge& e = edges[k];
++ if (v0 == e.v[1] && v1 == e.v[0])
++ {
++ add = false;
++ edges[k] = edges[--numEdges];
++ break;
++ }
++ }
++ if (add)
++ edges[numEdges++] = Edge::make(v0, v1);
++ }
++ mFaces[i] = mFaces[--mNumFaces];
++ }
++ }
++
++ if (numEdges == 0)
++ return FLT_MAX;
++
++ if (mNumFaces > MAX_FACES - numEdges)
++ return FLT_MAX;
++
++ for (PxU32 i = 0; i < numEdges; ++i)
++ {
++ const Edge& e = edges[i];
++ PxU8 i0 = e.v[0], i1 = e.v[1], i2 = PxU8(mNumVerts - 1);
++
++ PxVec4 plane;
++ if (!makePlane(i0, i1, i2, plane))
++ return FLT_MAX;
++
++ mFaces[mNumFaces++] = Face::make(plane, i0, i1, i2);
++ }
++
++ return proj;
++ }
++
++ PX_CUDA_CALLABLE
++ void computeClosest()
++ {
++ PX_ASSERT(mNumFaces > 0);
++
++ PxU32 closest = PxU32(-1);
++ PxReal closestW = -FLT_MAX;
++ for (PxU32 i = 0; i < mNumFaces; ++i)
++ {
++ const PxReal w = mFaces[i].plane.w;
++ if (w > closestW)
++ {
++ closest = i;
++ closestW = w;
++ }
++ }
++
++ PX_ASSERT(closest != PxU32(-1));
++
++ mClosest = mFaces[closest];
++
++ if (mNumFaces == 2)
++ mBest = mClosest;
++ }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxReal getDist() const
++ { return mClosest.plane.w; }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxVec3 getDir() const
++ { return mClosest.plane.getXYZ(); }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxReal getBestDist() const
++ { return mBest.plane.w; }
++
++ PX_CUDA_CALLABLE PX_INLINE
++ PxVec3 getBestDir() const
++ { return mBest.plane.getXYZ(); }
++
++ PX_CUDA_CALLABLE
++ void computePoints(PxVec3& pointA, PxVec3& pointB)
++ {
++ const Vert v0 = mVerts[mBest.v[0]], v1 = mVerts[mBest.v[1]], v2 = mVerts[mBest.v[2]];
++ const PxVec3 va0 = mConvexA.supportVertex(v0.aI), va1 = mConvexA.supportVertex(v1.aI), va2 = mConvexA.supportVertex(v2.aI);
++ const PxVec3 vb0 = mConvexB.supportVertex(v0.bI), vb1 = mConvexB.supportVertex(v1.bI), vb2 = mConvexB.supportVertex(v2.bI);
++ const PxVec3 a = va0 - vb0, b = va1 - vb1, c = va2 - vb2;
++ const PxVec3 abc = (b - a).cross(c - a);
++ const PxReal iabc = 1.0f / abc.magnitudeSquared();
++ const PxVec3 pabc = abc * abc.dot(a) * iabc;
++ const PxReal tbc = abc.dot((b - pabc).cross(c - pabc));
++ const PxReal tca = abc.dot((c - pabc).cross(a - pabc));
++ const PxReal tab = abc.dot((a - pabc).cross(b - pabc));
++ const PxReal sa = tbc * iabc, sb = tca * iabc, sc = tab * iabc;
++ pointA = va0 * sa + va1 * sb + va2 * sc;
++ pointB = vb0 * sa + vb1 * sb + vb2 * sc;
++ }
++
++ private:
++
++ static const PxU32 MAX_VERTS = 32;
++ static const PxU32 MAX_FACES = 2 * MAX_VERTS - 4;
++ static const PxU32 MAX_EDGES = MAX_VERTS + MAX_FACES - 2;
++
++ struct Vert
++ {
++ PxU8 aI, bI;
++
++ PX_CUDA_CALLABLE PX_INLINE
++ static Vert make(PxU8 aI, PxU8 bI)
++ { Vert v; v.aI = aI; v.bI = bI; return v; }
++ };
++
++ struct Face
++ {
++ PxVec4 plane;
++ PxU8 v[3];
++
++ PX_CUDA_CALLABLE PX_INLINE
++ static Face make(const PxVec4& plane, PxU8 v0, PxU8 v1, PxU8 v2)
++ { Face f; f.plane = plane; f.v[0] = v0; f.v[1] = v1; f.v[2] = v2; return f; }
++ };
++
++ struct Edge
++ {
++ PxU8 v[2];
++
++ PX_CUDA_CALLABLE PX_INLINE
++ static Edge make(PxU8 v0, PxU8 v1)
++ { Edge e; e.v[0] = v0; e.v[1] = v1; return e; }
++ };
++
++ Convex<Support>& mConvexA;
++ Convex<Support>& mConvexB;
++ Vert mVerts[MAX_VERTS];
++ Face mFaces[MAX_FACES];
++ PxU32 mNumVerts, mNumFaces;
++ Face mClosest, mBest;
++ PxReal mProj;
++ };
++
++ template <typename Support>
++ PX_CUDA_CALLABLE
++ static PxReal computeEpaDepth(const Support& a, const Support& b, const PxTransform& poseA, const PxTransform& poseB,
++ PxVec3& pointA, PxVec3& pointB, PxVec3& axis)
++ {
++ Convex<Support> convexA(a, poseA), convexB(b, poseB);
++ EpaDepth<Support> epa(convexA, convexB);
++ if (!epa.isValid())
++ return FLT_MAX;
++
++ const PxReal distEps = 0.01f * epa.getAccuracy();
++
++ while (true)
++ {
++ epa.computeClosest();
++
++ const PxReal closestDist = epa.getDist();
++ const PxVec3 dir = epa.getDir();
++
++ const PxReal proj = epa.addPoint(dir);
++
++ if (proj >= closestDist - distEps)
++ {
++ PxVec3 pA, pB;
++ epa.computePoints(pA, pB);
++ pointA = pA; pointB = pB;
++ axis = -epa.getBestDir();
++ return epa.getBestDist();
++ }
++ }
++ }
++ }
++ }
++ }
++
++ #endif