How to Calculate Surface Normals for Generated Geometry I have a class that generates a 3D shape based on inputs from the calling code. The inputs are things like length, depth, arc, etc. My code generates the geometry perfectly, however I am running into trouble when calculating the surface normals. When lit, my shape has very bizarre coloring/texture from the incorrect surface normals that are being calculated. From all my research I believe my math is correct, it seems that something is wrong with my technique or method.

At a high level how does one go about programmatically calculating the surface normals for a generated shape? I am using Swift/SceneKit on iOS for my code but a generic answer is fine.

I have two arrays that represent my shape. One is an array of 3D points that represents the vertices that make up the shape. The other array is a list of indexes of the first array that map the vertices into triangles. I need to take that data and generate a 3rd array that is a set of surface normals that aid in the lighting of the shape. (see SCNGeometrySourceSemanticNormal in SceneKit`)

The list of vertices and indexes is always different depending on the inputs to the class so I cannot pre-calculate or hard code the surface normals.
I see mainly three ways of computing normals for a generated shape.
Analytic normals
In some cases you have enough information about the surface to generate the normals. For example, the normal of any point on a sphere is trivial to compute. Put simply, when you know the derivative of the function, you also know the normal.
If your case is narrow enough to allow you to use analytic normals, they will probably give the best result in terms of precision. The technique doesn't scale too well though: if you also need to handle cases where you cannot use analytic normals, it may be easier to keep the technique that handles the general case and drop the analytic one altogether.
Vertex normals
The cross product of two vectors gives a vector perpendicular to the plane they belong to. So getting the normal of a triangle is straightforward:
vec3 computeNormal(vec3 a, vec3 b, vec3 c)
{
    return normalize(crossProduct(b - a, c - a));
}

Moreover, in the above example, the length of the cross product is proportional to the area inside abc. So the smoothed normal at a vertex shared by several triangles can be computed by summing up the cross products and normalizing as a last step, thus weighting each triangle by its area.
vec3 computeNormal(vertex a)
{
    vec3 sum = vec3(0, 0, 0);
    list&lt;vertex&gt; adjacentVertices = getAdjacentVertices(a);
    for (int i = 1; i &lt; adjacentVertices; ++i)
    {
        vec3 b = adjacentVertices[i - 1];
        vec3 c = adjacentVertices[i];
        sum += crossProduct(b - a, c - a);
    }
    if (norm(sum) == 0)
    {
        // Degenerate case
        return sum;
    }
    return normalize(sum);
}

If you are working with quads, there is a nice trick you can use: for a quad abcd, use crossProduct(c - a, d - b) and it will handle nicely cases where the quad is in fact a triangle.
Iñigo quilez wrote a few short articles on the topic: clever normalization of a mesh, and normal and area of n sided polygons.
Normals from partial derivatives
Normals can be computed in the fragment shader from the partial derivatives. The math behind is the same, except this time it is done in screen space. This article by Angelo Pesce describes the technique: Normals without normals.