I'm coding the method that inverts a matrix through Gaussian elimination. I've coded everything assuming that there are no zeros on the diagonal.
For the situation where a diagonal element is zero, what should I do? It's certain that some matrices with zero on the diagonal are still invertible, so how should I invert them by Gaussian elimination?
Create an "augmented matrix" of the form $[A|I]$ and use elementary row operations to put this into reduced row echelon form. The result will be $[I|A^{-1}]$ if $A$ is an invertible matrix.
Let's do an example, mainly for the benefit of future Readers. The OP mentions a "method [to] invert a matrix through Gaussian elimination ... assuming that there are no zero on the diagonal." It follows that the code already implements two of the three elementary row operations, multiplying a row by a nonzero scalar and adding (equiv. subtracting) a multiple of one row to (from) another row. What else is needed when we confront a zero "on the diagonal" (i.e. in a matrix position where we want a leading one) is the operation of exchanging two rows.
Using this operation does not "break" the matrix any more than do the other two kinds of elementary row operations. Yes, it changes the matrix, but the key here is that we want to change a matrix into its reduced row echelon form.
In particular the goal of inverting matrix $A$ will be achieved by finding the reduced row echelon form of the augmented matrix $[A|I]$. The required form can be found by an algorithmically determined sequence of elementary row operations, but the sequence is not unique. When no row exchanges are used, the method is called "Gaussian elimination without pivoting". This is possible precisely when the leading principal minors are nonzero, as this is another way of saying that zeros will not appear "on the diagonal," where we want the reduced row echelon form to have leading ones.
That is, elimination algorithms which use row exchanges commonly call those steps "pivoting". How hard one works to choose the "pivot" row is subject to a tradeoff in numeric stability, but in what follows we will present the case of partial pivoting (PDF). This amounts to taking the columns in whatever order they are presented (full pivoting meaning a strategy that might introduce leading ones in a wider consideration of all columns not yet reduced).
For our example we take:
$$ A = \begin{bmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{bmatrix} $$
Thus the augmented matrix we will reduce is $[A|I]$:
$$ \begin{bmatrix} \textbf 0 & 1 & 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 0 & 1 \end{bmatrix} $$
A reasonable rule when a zero appears where a leading one is wanted is to swap the row with the zero and one below it with an entry of largest maximum absolute value in the same column. Here we immediately start with a zero in the $(1,1)$ position, and below that both entries in that first column have the maximum absolute value $1$. So let's exchange the first and the second rows:
$$ \begin{bmatrix} \textbf 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 1 \end{bmatrix} $$
It is unnecessary to "scale" the leading entry in the current row (we already have a leading one). We proceed to clear the entries below it using addition/subtraction of a multiple of the first row to each other rows:
$$ \begin{bmatrix} \textbf 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 1 & 0 & 0 \\ 0 & 1 & -1 & 0 & -1 & 1 \end{bmatrix} $$
Having reduced the first column, we turn our attention to the second column. Fortuitously the $(2,2)$ entry is already a leading one in that row, and we proceed to add multiples of that second row to clear out the other entries in the second column:
$$ \begin{bmatrix} 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & \textbf 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & -2 & -1 & -1 & 1 \end{bmatrix} $$
Finally we scale the leading entry in the third row, so that the $(3,3)$ value becomes a leading one:
$$ \begin{bmatrix} 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & \textbf 1 & 1/2 & 1/2 & -1/2 \end{bmatrix} $$
Adding multiples of this new third row allows us to zero-out the other entries in the third column:
$$ \begin{bmatrix} 1 & 0 & 0 & -1/2 & 1/2 & 1/2 \\ 0 & 1 & 0 & 1/2 & -1/2 & 1/2 \\ 0 & 0 & \textbf 1 & 1/2 & 1/2 & -1/2 \end{bmatrix} $$
Now we have the reduced row echelon form, $[I|A^{-1}]$. To verify we can multiply:
$$ \begin{bmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{bmatrix} \; \begin{bmatrix} -1/2 & 1/2 & 1/2 \\ 1/2 & -1/2 & 1/2 \\ 1/2 & 1/2 & -1/2 \end{bmatrix} = I $$