Building indexed vertice/face structure.

I have a piece of code with which I am extracting triangles from a shape like so ... 

std::vector<double> points;   
// complete winding of points in triangle order.  
for (TopExp_Explorer exp (aShape, TopAbs_FACE); exp.More(); exp.Next())
  TriangleAccessor aTool (TopoDS::Face (exp.Current()));
  for (int iTri = 1; iTri <= aTool.NbTriangles(); iTri++)
    gp_Vec aNorm;
	gp_Pnt aPnt1, aPnt2, aPnt3;
    aTool.GetTriangle (iTri, aNorm, aPnt1, aPnt2, aPnt3);		
    points.push_back( aPnt1.X() ); points.push_back( aPnt1.Y() ); points.push_back( aPnt1.Z() );
    points.push_back( aPnt2.X() ); points.push_back( aPnt2.Y() ); points.push_back( aPnt2.Z() );
    points.push_back( aPnt3.X() ); points.push_back( aPnt3.Y() ); points.push_back( aPnt3.Z() ); 			

After this I am then building an indexed face list by checking those vertices to see which faces share them.

I see after a bit of coding that this is a naive approach because rounding issues will cause degenerate cases and the index information about each polygon is lost.

Now I am wondering if it is also possible using a topological explorer to extract edges that are shared from the shape and -then- relate those back to each triangle that I extract?

-OR- does a better mechanism exist for me to extract a properly indexed list of polygons and vertices?

What I am after is a mesh from this extraction that will test closed and manifold.

This is the code for the TriangleAccessor based on this ( ) 

// Auxiliary tools
  // Tool to get triangles from triangulation taking into account face
  // orientation and location
  class TriangleAccessor
    TriangleAccessor (const TopoDS_Face& aFace)
      TopLoc_Location aLoc;
      myPoly = BRep_Tool::Triangulation (aFace, aLoc);
      myTrsf = aLoc.Transformation();
      myNbTriangles = (myPoly.IsNull() ? 0 : myPoly->Triangles().Length());
      myInvert = (aFace.Orientation() == TopAbs_REVERSED);
      if (myTrsf.IsNegative())
        myInvert = ! myInvert;

    int NbTriangles () const { return myNbTriangles; } 

    // get i-th triangle and outward normal
    void GetTriangle (int iTri, gp_Vec &theNormal, gp_Pnt &thePnt1, gp_Pnt &thePnt2, gp_Pnt &thePnt3)
      // get positions of nodes
      int iNode1, iNode2, iNode3;
      myPoly->Triangles()(iTri).Get (iNode1, iNode2, iNode3); 
      thePnt1 = myPoly->Nodes()(iNode1);
      thePnt2 = myPoly->Nodes()(myInvert ? iNode3 : iNode2);
      thePnt3 = myPoly->Nodes()(myInvert ? iNode2 : iNode3);	

			int a = iNode1; 
			int b = (myInvert ? iNode3 : iNode2); 
			int c = (myInvert ? iNode2 : iNode3);

      // apply transormation if not identity
      if (myTrsf.Form() != gp_Identity)
        thePnt1.Transform (myTrsf);
        thePnt2.Transform (myTrsf);
        thePnt3.Transform (myTrsf);

      // calculate normal
      theNormal = (thePnt2.XYZ() - thePnt1.XYZ()) ^ (thePnt3.XYZ() - thePnt1.XYZ());
      Standard_Real aNorm = theNormal.Magnitude();
      if (aNorm > gp::Resolution())
        theNormal /= aNorm;

    Handle(Poly_Triangulation) myPoly;
    gp_Trsf myTrsf;
    int myNbTriangles;
    bool myInvert;


Would those nodes that are accessed during the triangle extraction also match nodes during an edge extraction. That is could I use the edge data to understand which edge on each triangle is shared with the edge on another triangle. 

Thanks !