Convex convex SAT

 procedure CollideConvexHullWithConvexHull(ShapeA,ShapeB:TEnginePhysicsShapeConvexHull);
 const kTolerance=0.005; // Skip near parallel edges: |Ea x Eb| = sin(alpha) * |Ea| * |Eb|
  function IsMinkowskiFace(const A,B,B_x_A,C,D,D_x_C:TVector3):boolean;
  var CBA,DBA,ADC,BDC:single;
  begin
   // Test if arcs AB and CD intersect on the unit sphere
   CBA:=Vector3Dot(C,B_x_A);
   DBA:=Vector3Dot(D,B_x_A);
   ADC:=Vector3Dot(A,D_x_C);
   BDC:=Vector3Dot(B,D_x_C);
   result:=((CBA*DBA<0.0)) and ((ADC*BDC)<0.0) and ((CBA*BDC)>0.0);
  end;
  function TestEarlyFaceDirection(const HullA,HullB:TEnginePhysicsShapeConvexHull;var FaceQuery:TEnginePhysicsContactFaceQuery):boolean;
  var Plane:TPlane;
      Transform:TMatrix4x4;
  begin
   Transform:=Matrix4x4TermMul(HullA.WorldTransform,HullB.InverseWorldTransform);
   Plane:=PlaneFastTransform(HullA.ConvexHull.Faces[FaceQuery.Index].Plane,Transform);
   FaceQuery.Separation:=PlaneVectorDistance(Plane,HullB.GetLocalSupport(Vector3Neg(Plane.Normal)));
   result:=FaceQuery.Separation>0.0;
  end;
  function TestEarlyEdgeDirection(const HullA,HullB:TEnginePhysicsShapeConvexHull;var EdgeQuery:TEnginePhysicsContactEdgeQuery):boolean;
  var EdgeA,EdgeB:PEnginePhysicsConvexHullEdge;
      L:single;
      CenterA,Pa,Qa,Ea,Ua,Va,Pb,Qb,Eb,Ub,Vb,Ea_x_Eb,Normal:TVector3;
      Transform:TMatrix4x4;
  begin
   result:=false;
   Transform:=Matrix4x4TermMul(HullA.WorldTransform,HullB.InverseWorldTransform);
   CenterA:=HullA.GetCenter(Transform);
   EdgeA:=@HullA.ConvexHull.Edges[EdgeQuery.IndexA];
   Pa:=Vector3TermMatrixMul(HullA.ConvexHull.Vertices[EdgeA^.Vertices[0]].Position,Transform);
   Qa:=Vector3TermMatrixMul(HullA.ConvexHull.Vertices[EdgeA^.Vertices[1]].Position,Transform);
   Ea:=Vector3Sub(Qa,Pa);
   Ua:=Vector3Norm(Vector3TermMatrixMulBasis(HullA.ConvexHull.Faces[EdgeA^.Faces[0]].Plane.Normal,Transform));
   Va:=Vector3Norm(Vector3TermMatrixMulBasis(HullA.ConvexHull.Faces[EdgeA^.Faces[1]].Plane.Normal,Transform));
   EdgeB:=@HullB.ConvexHull.Edges[EdgeQuery.IndexB];
   Pb:=HullB.ConvexHull.Vertices[EdgeB^.Vertices[0]].Position;
   Qb:=HullB.ConvexHull.Vertices[EdgeB^.Vertices[1]].Position;
   Eb:=Vector3Sub(Qb,Pb);
   Ub:=HullB.ConvexHull.Faces[EdgeB^.Faces[0]].Plane.Normal;
   Vb:=HullB.ConvexHull.Faces[EdgeB^.Faces[1]].Plane.Normal;
   if IsMinkowskiFace(Ua,Va,Vector3Neg(Ea),Vector3Neg(Ub),Vector3Neg(Vb),Vector3Neg(Eb)) then begin
    // Build search direction
    Ea_x_Eb:=Vector3Cross(Ea,Eb);

    // Skip near parallel edges: |Ea x Eb| = sin(alpha) * |Ea| * |Eb|
    L:=Vector3Length(Ea_x_Eb);
    if L<(sqrt(Vector3LengthSquared(Ea)*Vector3LengthSquared(Eb))*kTolerance) then begin
     result:=false;
     exit;
    end;

    // Assure consistent normal orientation (here: HullA -> HullB)
    Normal:=Vector3ScalarMul(Ea_x_Eb,1.0/L);
    if Vector3Dot(Normal,Vector3Sub(Pa,CenterA))<0.0 then begin
     Normal:=Vector3Neg(Normal);
    end;

    // s = Dot(Normal, Pb) - d = Dot(Normal, Pb) - Dot(Normal, Pa) = Dot(Normal, Pb - Pa)
    EdgeQuery.Separation:=Vector3Dot(Normal,Vector3Sub(Pb,Pa));
    if EdgeQuery.Separation>0.0 then begin
     result:=true;
    end;

   end;
  end;
  procedure QueryFaceDirections(const HullA,HullB:TEnginePhysicsShapeConvexHull;out OutFaceQuery:TEnginePhysicsContactFaceQuery);
  var MaxIndex,Index:longint;
      MaxSeparation,Separation:single;
      Plane:TPlane;
      Transform:TMatrix4x4;
  begin
   Transform:=Matrix4x4TermMul(HullA.WorldTransform,HullB.InverseWorldTransform);
   MaxIndex:=-1;
   MaxSeparation:=-3.4e+38;
   for Index:=0 to HullA.ConvexHull.CountFaces-1 do begin
    Plane:=PlaneFastTransform(HullA.ConvexHull.Faces[Index].Plane,Transform);
    Separation:=PlaneVectorDistance(Plane,HullB.GetLocalSupport(Vector3Neg(Plane.Normal)));
    if (Index=0) or (MaxSeparation<Separation) then begin
     MaxSeparation:=Separation;
     MaxIndex:=Index;
     if MaxSeparation>0.0 then begin
      break;
     end;
    end;
   end;
   OutFaceQuery.Index:=MaxIndex;
   OutFaceQuery.Separation:=MaxSeparation;
  end;
  procedure QueryEdgeDirections(const HullA,HullB:TEnginePhysicsShapeConvexHull;out OutEdgeQuery:TEnginePhysicsContactEdgeQuery);
  var EdgeA,EdgeB:PEnginePhysicsConvexHullEdge;
      IndexA,IndexB,MaxIndexA,MaxIndexB:longint;
      MaxSeparation,Separation,L:single;
      CenterA,Pa,Qa,Ea,Ua,Va,Pb,Qb,Eb,Ub,Vb,Ea_x_Eb,Normal:TVector3;
      Transform:TMatrix4x4;
      First:boolean;
  begin
   MaxIndexA:=-1;
   MaxIndexB:=-1;
   MaxSeparation:=-3.4e+38;
   Transform:=Matrix4x4TermMul(HullA.WorldTransform,HullB.InverseWorldTransform);
   CenterA:=HullA.GetCenter(Transform);
   First:=true;
   for IndexA:=0 to HullA.ConvexHull.CountEdges-1 do begin
    EdgeA:=@HullA.ConvexHull.Edges[IndexA];
    Pa:=Vector3TermMatrixMul(HullA.ConvexHull.Vertices[EdgeA^.Vertices[0]].Position,Transform);
    Qa:=Vector3TermMatrixMul(HullA.ConvexHull.Vertices[EdgeA^.Vertices[1]].Position,Transform);
    Ea:=Vector3Sub(Qa,Pa);
    Ua:=Vector3Norm(Vector3TermMatrixMulBasis(HullA.ConvexHull.Faces[EdgeA^.Faces[0]].Plane.Normal,Transform));
    Va:=Vector3Norm(Vector3TermMatrixMulBasis(HullA.ConvexHull.Faces[EdgeA^.Faces[1]].Plane.Normal,Transform));
    for IndexB:=0 to HullB.ConvexHull.CountEdges-1 do begin
     EdgeB:=@HullB.ConvexHull.Edges[IndexB];
     Pb:=HullB.ConvexHull.Vertices[EdgeB^.Vertices[0]].Position;
     Qb:=HullB.ConvexHull.Vertices[EdgeB^.Vertices[1]].Position;
     Eb:=Vector3Sub(Qb,Pb);
     Ub:=HullB.ConvexHull.Faces[EdgeB^.Faces[0]].Plane.Normal;
     Vb:=HullB.ConvexHull.Faces[EdgeB^.Faces[1]].Plane.Normal;
     if IsMinkowskiFace(Ua,Va,Vector3Neg(Ea),Vector3Neg(Ub),Vector3Neg(Vb),Vector3Neg(Eb)) then begin
      // Build search direction
      Ea_x_Eb:=Vector3Cross(Ea,Eb);

      // Skip near parallel edges: |Ea x Eb| = sin(alpha) * |Ea| * |Eb|
      L:=Vector3Length(Ea_x_Eb);
      if L<(sqrt(Vector3LengthSquared(Ea)*Vector3LengthSquared(Eb))*kTolerance) then begin
       continue;
      end;

      // Assure consistent normal orientation (here: HullA -> HullB)
      Normal:=Vector3ScalarMul(Ea_x_Eb,1.0/L);
      if Vector3Dot(Normal,Vector3Sub(Pa,CenterA))<0.0 then begin
       Normal:=Vector3Neg(Normal);
      end;

      // s = Dot(Normal, Pb) - d = Dot(Normal, Pb) - Dot(Normal, Pa) = Dot(Normal, Pb - Pa)
      Separation:=Vector3Dot(Normal,Vector3Sub(Pb,Pa));
      if First or (MaxSeparation<Separation) then begin
       First:=false;
       MaxSeparation:=Separation;
       MaxIndexA:=IndexA;
       MaxIndexB:=IndexB;
       if MaxSeparation>0.0 then begin
        break;
       end;
      end;

     end;
    end;
   end;
   OutEdgeQuery.IndexA:=MaxIndexA;
   OutEdgeQuery.IndexB:=MaxIndexB;
   OutEdgeQuery.Separation:=MaxSeparation;
  end;
  function GetEdgeContact(var CA,CB,Normal:TVector3;const PA,QA,PB,QB,DirectionAB:TVector3):boolean;
  var DA,DB,r:TVector3;
      a,e,f,c,b,d,TA,TB:single;
  begin
   DA:=Vector3Sub(QA,PA);
   DB:=Vector3Sub(QB,PB);
   r:=Vector3Sub(PA,PB);
     a:=Vector3LengthSquared(DA);
     e:=Vector3LengthSquared(DB);
     f:=Vector3Dot(DB,r);
     c:=Vector3Dot(DA,r);
   b:=Vector3Dot(DA,DB);
   d:=(a*e)-sqr(b);
   if (abs(d)>EPSILON) and (abs(e)>EPSILON) then begin

    TA:=((b*f)-(c*e))/d;
    TB:=((b*TA)+f)/e;

    CA:=Vector3Add(PA,Vector3ScalarMul(DA,TA));
    CB:=Vector3Add(PB,Vector3ScalarMul(DB,TB));

    Normal:=Vector3Norm(Vector3Cross(DA,DB));

    // Assure consistent normal orientation (here: HullA -> HullB)
    if Vector3Dot(Normal,DirectionAB)<0.0 then begin
     Normal:=Vector3Neg(Normal);
    end;

    result:=true;

   end else begin
    result:=false;
   end;
  end;
  function FindIncidentFaceIndex(const ReferenceHull:TEnginePhysicsShapeConvexHull;const ReferenceFaceIndex:longint;const IncidentHull:TEnginePhysicsShapeConvexHull):longint;
  var i:longint;
      MinDot,Dot:single;
      ReferenceNormal:TVector3;
  begin
   ReferenceNormal:=Vector3TermMatrixMulBasis(Vector3TermMatrixMulBasis(ReferenceHull.ConvexHull.Faces[ReferenceFaceIndex].Plane.Normal,
                                                                        ReferenceHull.WorldTransform),
                                              IncidentHull.InverseWorldTransform);
   result:=-1;
   MinDot:=3.4e+38;
   for i:=0 to IncidentHull.ConvexHull.CountFaces-1 do begin
    Dot:=Vector3Dot(ReferenceNormal,IncidentHull.ConvexHull.Faces[i].Plane.Normal);
        if MinDot>Dot then begin
         MinDot:=Dot;
     result:=i;
    end;
   end;
  end;
  procedure ClipFaceContactPoints(const ReferenceHull:TEnginePhysicsShapeConvexHull;const ReferenceFaceIndex:longint;const IncidentHull:TEnginePhysicsShapeConvexHull;const IncidentFaceIndex:longint);
  var Contact:PEnginePhysicsContact;
      ReferenceVertexIndex,OtherReferenceVertexIndex,IncidentVertexIndex,ClipVertexIndex:longint;
      ReferenceFace,IncidentFace:PEnginePhysicsConvexHullFace;
      ClipVertex,FirstClipVertex,EndClipVertex:PVector3;
      StartDistance,EndDistance,Distance:single;
      ClipVertices:array[0..2] of TEnginePhysicsConvexHullVertexList;
      ReferencePoint:TVector3;
      ReferenceWorldPlane,ReferenceEdgePlane:TPlane;
  begin

   ContactManager.CountTemporaryContacts[ThreadIndex]:=0;

   ReferenceFace:=@ReferenceHull.ConvexHull.Faces[ReferenceFaceIndex];
   ReferenceWorldPlane:=PlaneFastTransform(ReferenceFace^.Plane,ReferenceHull.WorldTransform);

   IncidentFace:=@IncidentHull.ConvexHull.Faces[IncidentFaceIndex];

   ClipVertices[0]:=ContactManager.ConvexHullVertexLists[ThreadIndex,0];
   ClipVertices[0].Clear;

   for IncidentVertexIndex:=0 to IncidentFace^.CountVertices-1 do begin
    ClipVertices[0].Add(Vector3TermMatrixMul(IncidentHull.ConvexHull.Vertices[IncidentFace^.Vertices[IncidentVertexIndex]].Position,IncidentHull.WorldTransform));
   end;

   ClipVertices[1]:=ContactManager.ConvexHullVertexLists[ThreadIndex,1];
   ClipVertices[1].Clear;

   OtherReferenceVertexIndex:=ReferenceFace^.CountVertices-1;
   for ReferenceVertexIndex:=0 to ReferenceFace^.CountVertices-1 do begin
    if ClipVertices[0].Count>=2 then begin
     ReferencePoint:=Vector3TermMatrixMul(ReferenceHull.ConvexHull.Vertices[ReferenceFace^.Vertices[ReferenceVertexIndex]].Position,ReferenceHull.WorldTransform);
     ReferenceEdgePlane.Normal:=Vector3Neg(Vector3Norm(Vector3Cross(ReferenceWorldPlane.Normal,Vector3Sub(ReferencePoint,Vector3TermMatrixMul(ReferenceHull.ConvexHull.Vertices[ReferenceFace^.Vertices[OtherReferenceVertexIndex]].Position,ReferenceHull.WorldTransform)))));
     ReferenceEdgePlane.Distance:=-Vector3Dot(ReferenceEdgePlane.Normal,ReferencePoint);
     FirstClipVertex:=@ClipVertices[0].Vertices[ClipVertices[0].Count-1];
     EndClipVertex:=@ClipVertices[0].Vertices[0];
     StartDistance:=PlaneVectorDistance(ReferenceEdgePlane,FirstClipVertex^);
     for ClipVertexIndex:=0 to ClipVertices[0].Count-1 do begin
      EndClipVertex:=@ClipVertices[0].Vertices[ClipVertexIndex];
      EndDistance:=PlaneVectorDistance(ReferenceEdgePlane,EndClipVertex^);
      if StartDistance<0.0 then begin
       if EndDistance<0.0 then begin
        ClipVertices[1].Add(EndClipVertex^);
       end else begin
        ClipVertices[1].Add(Vector3Lerp(FirstClipVertex^,EndClipVertex^,StartDistance/(StartDistance-EndDistance)));
       end;
      end else if EndDistance<0.0 then begin
       ClipVertices[1].Add(Vector3Lerp(FirstClipVertex^,EndClipVertex^,StartDistance/(StartDistance-EndDistance)));
       ClipVertices[1].Add(EndClipVertex^);
      end;
      FirstClipVertex:=EndClipVertex;
      StartDistance:=EndDistance;
     end;
    end;
    if ClipVertices[1].Count=0 then begin
     exit;
    end else begin
     ClipVertices[2]:=ClipVertices[0];
     ClipVertices[0]:=ClipVertices[1];
     ClipVertices[1]:=ClipVertices[2];
     ClipVertices[1].Clear;
     OtherReferenceVertexIndex:=ReferenceVertexIndex;
    end;
   end;

   for ClipVertexIndex:=0 to ClipVertices[0].Count-1 do begin
    ClipVertex:=@ClipVertices[0].Vertices[ClipVertexIndex];
    Distance:=PlaneVectorDistance(ReferenceWorldPlane,ClipVertex^);
    if Distance<0.0 then begin
     if ContactManager.CountTemporaryContacts[ThreadIndex]<MAX_TEMPORARY_CONTACTS then begin
      Contact:=@ContactManager.TemporaryContacts[ThreadIndex,ContactManager.CountTemporaryContacts[ThreadIndex]];
      inc(ContactManager.CountTemporaryContacts[ThreadIndex]);
      Contact^.WorldPositions[0]:=ClipVertex^;
      Contact^.WorldPositions[1]:=ClipVertex^;
      Contact^.LocalPositions[0]:=Vector3TermMatrixMul(ClipVertex^,ShapeA.InverseWorldTransform);
      Contact^.LocalPositions[1]:=Vector3TermMatrixMul(ClipVertex^,ShapeB.InverseWorldTransform);
      Contact^.Penetration:=Distance;
      Contact^.FeaturePair.Key:=$ffffffff; // $ffffffff => nearest current-frame=>last-frame contact point search for warm starting
     end else begin
      break;
     end;
    end;
   end;

  end;
 var Contact:PEnginePhysicsContact;
     Iteration,ReferenceFaceIndex,IncidentFaceIndex:longint;
     EdgeA,EdgeB:PEnginePhysicsConvexHullEdge;
     PenetrationDepth:single;
     pa,pb,Normal:TVector3;
     ConvexShapes:array[0..1] of TEnginePhysicsConvexShape;
 begin

  Manifold.Persistent:=false;
  Manifold.CountContacts:=0;

  Iteration:=0;
  while Iteration<2 do begin

   ContactManager.CountTemporaryContacts[ThreadIndex]:=0;

   if Iteration=0 then begin

    if Manifold.HaveData then begin

     if (Manifold.FaceQueryAB.Index>=0) and (Manifold.FaceQueryAB.Separation>0.0) then begin
      if TestEarlyFaceDirection(ShapeA,ShapeB,Manifold.FaceQueryAB) then begin
       // Still existent seperating axis from last frame found, so exit!
       exit;
      end else begin
       // Reject the try to rebuild the contact manifold from last frame, and process a new full seperating axis test
       Iteration:=1;
       continue;
      end;
     end else if (Manifold.FaceQueryBA.Index>=0) and (Manifold.FaceQueryBA.Separation>0.0) then begin
      if TestEarlyFaceDirection(ShapeB,ShapeA,Manifold.FaceQueryBA) then begin
       // Still existent seperating axis from last frame found, so exit!
       exit;
      end else begin
       // Reject the try to rebuild the contact manifold from last frame, and process a new full seperating axis test
       Iteration:=1;
       continue;
      end;
     end else if ((Manifold.EdgeQuery.IndexA>=0) and (Manifold.EdgeQuery.IndexB>=0)) and (Manifold.EdgeQuery.Separation>0.0) then begin
      if TestEarlyEdgeDirection(ShapeA,ShapeB,Manifold.EdgeQuery) then begin
       // Still existent seperating axis from last frame found, so exit!
       exit;
      end else begin
       // Reject the try to rebuild the contact manifold from last frame, and process a new full seperating axis test
       Iteration:=1;
       continue;
      end;
     end else if ((Manifold.FaceQueryAB.Index<0) or (Manifold.FaceQueryAB.Separation>0.0)) or
                 ((Manifold.FaceQueryBA.Index<0) or (Manifold.FaceQueryBA.Separation>0.0)) or
                 (((Manifold.EdgeQuery.IndexA<0) or (Manifold.EdgeQuery.IndexB<0)) or (Manifold.EdgeQuery.Separation>0.0)) then begin
      // Reject the try to rebuild the contact manifold from last frame, and process a new full seperating axis test
      Iteration:=1;
      continue;
     end else begin
      if TestEarlyFaceDirection(ShapeA,ShapeB,Manifold.FaceQueryAB) then begin
       // Still existent seperating axis from last frame found, so exit!
       exit;
      end else if TestEarlyFaceDirection(ShapeB,ShapeA,Manifold.FaceQueryBA) then begin
       // Still existent seperating axis from last frame found, so exit!
       exit;
      end else if TestEarlyEdgeDirection(ShapeA,ShapeB,Manifold.EdgeQuery) then begin
       // Still existent seperating axis from last frame found, so exit!
       exit;
      end else if ((Manifold.EdgeQuery.IndexA>=0) and (Manifold.EdgeQuery.IndexB>=0)) and
                  ((Manifold.EdgeQuery.Separation>(Manifold.FaceQueryAB.Separation+kTolerance)) and (Manifold.EdgeQuery.Separation>(Manifold.FaceQueryBA.Separation+kTolerance))) then begin
       // Reject the try to rebuild the contact manifold from last frame, and process a new full seperating axis test
       Iteration:=1;
       continue;
      end else begin
       // Okay in this case, we can try to rebuild the contact manifold from last frame
      end;
     end;

    end else begin

     // We must process a full seperating axis test, since there are no last frame contact manifold data yet
     Iteration:=1;
     continue;

    end;

   end else begin

    Manifold.FaceQueryAB.Index:=-1;
    Manifold.FaceQueryAB.Separation:=3.4e+38;

    Manifold.FaceQueryBA.Index:=-1;
    Manifold.FaceQueryBA.Separation:=3.4e+38;

    Manifold.EdgeQuery.IndexA:=-1;
    Manifold.EdgeQuery.IndexB:=-1;
    Manifold.EdgeQuery.Separation:=3.4e+38;

    Manifold.HaveData:=true;

    QueryFaceDirections(ShapeA,ShapeB,Manifold.FaceQueryAB);
    if Manifold.FaceQueryAB.Separation>0.0 then begin
     exit;
    end;

    QueryFaceDirections(ShapeB,ShapeA,Manifold.FaceQueryBA);
    if Manifold.FaceQueryBA.Separation>0.0 then begin
     exit;
    end;

    QueryEdgeDirections(ShapeA,ShapeB,Manifold.EdgeQuery);
    if Manifold.EdgeQuery.Separation>0.0 then begin
     exit;
    end;

   end;

   if ((Manifold.EdgeQuery.IndexA>=0) and (Manifold.EdgeQuery.IndexB>=0)) and
      ((Manifold.EdgeQuery.Separation>(Manifold.FaceQueryAB.Separation+kTolerance)) and (Manifold.EdgeQuery.Separation>(Manifold.FaceQueryBA.Separation+kTolerance))) then begin

    // Edge contact

    Manifold.HaveData:=false;

    EdgeA:=@ShapeA.ConvexHull.Edges[Manifold.EdgeQuery.IndexA];
    EdgeB:=@ShapeB.ConvexHull.Edges[Manifold.EdgeQuery.IndexB];

    if GetEdgeContact(pa,
                      pb,
                      Normal,
                      Vector3TermMatrixMul(ShapeA.ConvexHull.Vertices[EdgeA^.Vertices[0]].Position,ShapeA.WorldTransform),
                      Vector3TermMatrixMul(ShapeA.ConvexHull.Vertices[EdgeA^.Vertices[1]].Position,ShapeA.WorldTransform),
                      Vector3TermMatrixMul(ShapeB.ConvexHull.Vertices[EdgeB^.Vertices[0]].Position,ShapeB.WorldTransform),
                      Vector3TermMatrixMul(ShapeB.ConvexHull.Vertices[EdgeB^.Vertices[1]].Position,ShapeB.WorldTransform),
                      Vector3Sub(ShapeA.GetCenter(ShapeA.WorldTransform),ShapeB.GetCenter(ShapeB.WorldTransform))) then begin
     PenetrationDepth:=Vector3Dot(Vector3Sub(pb,pa),Normal);
     if PenetrationDepth<0.0 then begin
      Manifold.Normal:=Vector3Neg(Normal);
      Contact:=@ContactManager.TemporaryContacts[ThreadIndex,ContactManager.CountTemporaryContacts[ThreadIndex]];
      inc(ContactManager.CountTemporaryContacts[ThreadIndex]);
      Contact^.WorldPositions[0]:=Vector3Avg(pa,pb);
      Contact^.WorldPositions[1]:=Contact^.WorldPositions[0];
      Contact^.LocalPositions[0]:=Vector3TermMatrixMul(Contact^.WorldPositions[0],ShapeA.InverseWorldTransform);
      Contact^.LocalPositions[1]:=Vector3TermMatrixMul(Contact^.WorldPositions[1],ShapeB.InverseWorldTransform);
      Contact^.Penetration:=PenetrationDepth;
      Contact^.FeaturePair.Key:=$ffffffff; // $ffffffff => nearest current-frame=>last-frame contact point search for warm starting
     end;
    end;

    if (Iteration>0) and (ContactManager.CountTemporaryContacts[ThreadIndex]=0) then begin
     // MPR to the rescue! This case should never happen, but secure is secure :-)
     ConvexShapes[0].Data:=ShapeA;
     ConvexShapes[0].Transform:=@ShapeA.WorldTransform;
     ConvexShapes[0].CenterFunction:=@ConvexShapeCenterFunctionShape;
     ConvexShapes[0].SupportFunction:=@ConvexShapeSupportFunctionShape;
     ConvexShapes[1].Data:=ShapeB;
     ConvexShapes[1].Transform:=@ShapeB.WorldTransform;
     ConvexShapes[1].CenterFunction:=@ConvexShapeCenterFunctionShape;
     ConvexShapes[1].SupportFunction:=@ConvexShapeSupportFunctionShape;
     if MPRPenetration(@ConvexShapes[0],@ConvexShapes[1],pa,pb,Normal,PenetrationDepth) then begin
      Manifold.Normal:=Vector3Neg(Normal);
      Contact:=@ContactManager.TemporaryContacts[ThreadIndex,ContactManager.CountTemporaryContacts[ThreadIndex]];
      inc(ContactManager.CountTemporaryContacts[ThreadIndex]);
      Contact^.WorldPositions[0]:=Vector3Avg(pa,pb);
      Contact^.WorldPositions[1]:=Contact^.WorldPositions[0];
      Contact^.LocalPositions[0]:=Vector3TermMatrixMul(Contact^.WorldPositions[0],ShapeA.InverseWorldTransform);
      Contact^.LocalPositions[1]:=Vector3TermMatrixMul(Contact^.WorldPositions[1],ShapeB.InverseWorldTransform);
      Contact^.Penetration:=-PenetrationDepth;
      Contact^.FeaturePair.Key:=$ffffffff; // $ffffffff => nearest current-frame=>last-frame contact point search for warm starting
     end;
    end;

   end else begin

    // Face contact

    // TODO: Is it the (or also a) correct way, that the reference face for
    // clipping in this implementation is always on HullA, and the clipping
    // incident face always on HullB? At least, in this case the contact
    // points are always directly automatically coherence and consistent
    // without plane reprojection (among other things for warm starting for
    // the nearest new->old feature-ID-key-free contact point search between
    // last and current frame). It seems, that Bullet does also in this way.
    // And this way implementation works pretty well so far. At least, I saw
    // no failure cases at my tests so far.
    //
    // The alternative solution should be that a face from HullB could also
    // the clip reference face, where the contact points would be projected
    // to the HullA in-this-case-now incident face then, but will be it the
    // more correct way? It seems, that the most other physics engines
    // (but not Bullet) does in this way, but these uses often contact pair
    // feature ID keys then.
    //
    // Or in other words, what are the pros&cons of the two possible
    // face-face-clipping implementation ways?

    if (Manifold.FaceQueryAB.Separation+EPSILON)>Manifold.FaceQueryBA.Separation then begin

     ReferenceFaceIndex:=Manifold.FaceQueryAB.Index;
     IncidentFaceIndex:=FindIncidentFaceIndex(ShapeA,Manifold.FaceQueryAB.Index,ShapeB);

     Manifold.Normal:=Vector3TermMatrixMulBasis(ShapeA.ConvexHull.Faces[Manifold.FaceQueryAB.Index].Plane.Normal,ShapeA.WorldTransform);

    end else begin

     ReferenceFaceIndex:=FindIncidentFaceIndex(ShapeB,Manifold.FaceQueryBA.Index,ShapeA);
     IncidentFaceIndex:=Manifold.FaceQueryBA.Index;

     // The HullB face normal must be corrected (negated) for a still consistent normal orientation
     Manifold.Normal:=Vector3Neg(Vector3TermMatrixMulBasis(ShapeB.ConvexHull.Faces[Manifold.FaceQueryBA.Index].Plane.Normal,ShapeB.WorldTransform));

    end;
    ClipFaceContactPoints(ShapeA,ReferenceFaceIndex,ShapeB,IncidentFaceIndex);

   end;

   if ContactManager.CountTemporaryContacts[ThreadIndex]>0 then begin
    // Contacts found, reduce these down to four contacts with the largest area
    Manifold.CountContacts:=ContactManager.ReduceContacts(pointer(@ContactManager.TemporaryContacts[ThreadIndex,0]),ContactManager.CountTemporaryContacts[ThreadIndex],pointer(@Manifold.Contacts[0]));
    exit;
   end else begin
    if Iteration=0 then begin
     // We must process a new full seperating axis test, since the last frame contact manifold could not rebuilt.
     inc(Iteration);
    end else begin
     // No contacts found
     exit;
    end;
   end;

  end;

 end;