For a differentiable surface, we can use calculus to find the normals of points on the surface. Unfortunately, a triangle mesh is not

differentiable. The technique that is generally applied to triangle meshes is called vertex normal averaging. The vertex normal n for

an arbitrary vertex v in a mesh is found by averaging the face normals of every polygon in the mesh that shares the vertex v. For

example, in Figure 1, four polygons in the mesh share the vertex v; thus, the vertex normal for v is given by:

Figure 1

The middle vertex is shared by the neighboring four polygons, so we approximate the middle vertex normal by averaging the four polygon

face normals.

<-- TO FIND Navg

Im familiar with finding face normal in a triangle with this formula,

But how do I find n0,n1,n2,n3 if there not in a triangle, but in the middle of a square as in Figure 1

The normal is the same for the entire plane where the 3 vertex positions (triangle) or 4 vertex positions (quad) lay in. So you do not compute the normal in the middle of a square, you do compute the normal for the said entire plane.

With the assumption that the quad given by 4 corners in order p0, p1, p2, p3, is planar and not degenerate, the difference vectors p1-p0 and p2-p0 can be used as u and v. Often also the "cross over" difference vectors p3-p1 and p2-p0 are used if a quad is given.

BTW: Computing a vertex normal is usually done not using the normals of the adjacent faces directly but using them weighted. The weights are often chosen w.r.t. the angle the face has at the vertex, so that smaller angles will have less effect.

In addition to haegarr's excellent post (thumbs up), depending on the effect to be achieved (you don't mention the purpose you intend), weighting of the adjacent face normals can also be proportional to the triangle/quad area.