I know this question is already answered, but here's a similar single-face based looping storage which is implemented in the OpenFOAM C++ library:
Each cell has an index (ID) in a cellList.
Two lists are defined for all faces: "face internal owner" and "face neighbour".
The length of both face lists corresponds to the number of internal faces in the mesh.
A face owner will be the cell with the lower ID in the cellList (opposite for face neighbour).
Boundary faces are written last, and they have outward oriented normals (from the solution domain), and of course, only one owner cell.
The face area normal is oriented so that it looks outwards from the owner cell to the neighbour cell.
This works well for e.g. flux calculation. The flux is evaluated once per face, and it is added to the sum of total faces for the owner cells, and deducted from the neighbour cells (the summation/deduction is decided based on the orientation of the face area normal). The boundary faces are sorted and stored at the bottom of the face list, allowing boundary conditions to be defined as slices of the face list (begining label, end label of the boundary patch),simplifying thus the implementation of the boundary conditions, as well as inhancing efficiency of the updating process for the boundary conditions, since it is relying on the solution provided by the operations on internal faces.
Since the boundary faces are agglomerated into patches, the inter-process communication is defined for coupled (processor) patches, and pre-defined. This means that as soon as there is a loop over the boundary mesh, the top level access functions envoke wrapped MPI calls, making such code "automatically" parallelized, if it relies on the above explained face-based connectivity.