I'll add this to the other answers, which do a very good job of explaining why this very irregular shape is not inconsistent with the Moon's original much "rounder" shape.
This could be done with some geometry and math, but I was lazy and did a rendering in Blender. I made a spheroid and added some bumps, then added a Sun lamp above and screen far below. I can adjust the half-angle of the sun lamp. This is not intended to be an exact simulation, but once you see the GIF the principle described in words in the other two answers is reinforced.
Note the appearance of sharp corners on the shape of the "umbra" that do not exist on the original shape!
In order to delineate the "umbra" where the light reaches zero, I just turned up the brightness of the Sun to 100
.
below: The original shape, moved down into the plane of the screen so that the cross-section shape is clear.
below: Selected frames, showing a nearly square shaped "umbra" with the appearance of sharp corners that do not exist on the original shape!
Here's the Python script I used in Blender to make the shape. Animation is easily scripted also, I just did it manually this time because it was quicker.
import bpy
import numpy as np
halfpi, pi, twopi, fourpi = [f*np.pi for f in [0.5, 1, 2, 4]]
half_nth = 50
nth, nph = 2*half_nth + 1, 200
th = np.linspace(-halfpi, halfpi, nth+2)[1:-1]
phi = np.linspace(0, twopi, nph+1)[:-1]
n1, n2 = 8, 32
thing = 0.5*(1.+np.cos(np.linspace(0, pi, n2+1-n1)))[:,None]
ramp = np.zeros_like(th)[:,None]
ramp[half_nth+n1 : half_nth+(n2+1)] = thing
ramp[half_nth-n2 : half_nth-n1+1 ] = thing[::-1]
ramp[half_nth-n1+1 : half_nth+n1 ] = 1.0
nwiggle = 4
awiggle = 0.05
wiggle = 1. + awiggle * np.cos(nwiggle*phi) * ramp
X = wiggle * np.cos(phi)[None,:] * np.cos(th)[:,None]
Y = wiggle * np.sin(phi)[None,:] * np.cos(th)[:,None]
Z = np.ones_like(phi)[None,:] * np.sin(th)[:,None]
verts = [tuple(thing) for thing in zip(X.flatten(), Y.flatten(), Z.flatten())]
faces = []
for ith in range(nth-1):
for iph in range(nph):
v1 = iph + nph*ith
v2 = (iph+1)%nph + nph*ith
v3 = (iph+1)%nph + nph*(ith+1)
v4 = iph + nph*(ith+1)
faces.append((v1, v2, v3, v4))
bot = [tuple([iph for iph in range(nph)])]
top = [tuple([iph + nph*(nth-1) for iph in range(nph)])]
faces += bot + top
me = bpy.data.meshes.new('wow')
ob = bpy.data.objects.new('wow', me)
ob.location = (0.0, 0.0, 102.0)
bpy.context.scene.objects.link(ob)
me.from_pydata(verts,[],faces)
bpy.data.objects['wow'].select = False
bpy.data.objects['wow'].select = True
bpy.ops.object.shade_smooth()