Be aware that if you will access some surface<Type>Field given the faces, you need to check whether or not the face label is internal. If not, then you need to access the value of this boundary face through the boundary field.

basically a need to smooth alpha on every cell and i would like to do that by using alpha on the face times area divided by the total area. So i need to sum alpha on the 4 faces. The divergence is a difference between values on the faces divided by deltax and i dont know if is the same that i need.

I decided to post my question inside this thread, as it is also somehow connected with looping over faces in a given cell. What I try to achieve is to get the face normal vector at every face in the cell, regardless it is a boundary cell or internal.
I produced the following code...:

...but the values of normal vectors I received are all of positive sign, all the components are of positive sign (on hexahedral mesh). I think in three normal vectors there should be at least one component of negative sign... Do you have any ideas?

Yes, I solved the problem. The solution I got was correct. I.e. the components of normal vectors for an arbitrary cell can be in general of the same sign.

The reason for this are the normal vector direction rules. In OpenFoam mesh normal vectors:
- at boundary faces point out of the domain
- at internal faces they point from the cell of lower global ID number to higher.