When they solve $Av=-b$ the assumption is that $A$ is non-singular. Whenever $A$ is singular the problem can either have no solution, or infinitely many solutions. In their case I believe it is the latter. Remember that $v = (p,s)$. What they say is that they set the solution for $p$ to be the edge midpoint: $p=\frac{p_1+p_2}{2}$ and then plug that into the equation to reduce the system:
\begin{align}
-b = Av \iff -\begin{bmatrix} b_1\\ b_2\end{bmatrix} =
\begin{bmatrix}
A_{11} & A_{12}\\
A_{21} & I
\end{bmatrix}
\begin{bmatrix}
p \\ s \end{bmatrix} =
\begin{bmatrix} A_{11} \\ A_{21}\end{bmatrix}p+ \begin{bmatrix} A_{12} \\ I \end{bmatrix} s.
\end{align}
If you assume the system is consistent, then $s = -b_2 - A_{21}p$. If the system is inconsistent, then you can try to find the solution that least violates it as:
$$(A_{12}^TA_{12}+I)s = -A_{12}^T(b_1+A_{11}p) - (b_2+A_{21}p).$$
I suppose the system ought to be consistent though.
Edit:
Here's a more general approach as to how you can factorize such problems if you assume $A_{22}$ is an invertible matrix (it is essentially block Gaussian elimination):
\begin{align}
\begin{bmatrix} b_1 \\ b_2\end{bmatrix} &=
\begin{bmatrix}
A_{11} & A_{12} \\ A_{21} & A_{22}
\end{bmatrix}
\begin{bmatrix}
x_1\\ x_2
\end{bmatrix}
\\
\begin{bmatrix}
b_1 \\ b_2
\end{bmatrix}
&= \begin{bmatrix}
A_{11}x_1 + A_{12}x_2 \\
A_{21}x_1 + A_{22}x_2
\end{bmatrix}
\\
\begin{bmatrix}
b_1\\
A_{22}^{-1}(b_2 - A_{21}x_1)
\end{bmatrix}
&=
\begin{bmatrix}
A_{11}x_1 + A_{12}x_2
\\
x_2
\end{bmatrix}
\\
\begin{bmatrix}
b_1 - A_{12}A_{22}^{-1}b_2 \\
A_{22}^{-1}(b_2 - A_{21}x_1)
\end{bmatrix}
&=
\begin{bmatrix}
(A_{11}-A_{12}A_{22}^{-1}A_{21})x_1 \\ x_2
\end{bmatrix}
\end{align}
In your case $A_{22} = I$ so the inverse is just $I$. If the Schur complement $[A/A_{22}]= A_{11}-A_{12}A_{22}^{-1}A_{21}$ is invertible, then
$$x_1 = [A/A_{22}]^{-1}(b_1 - A_{12}A_{22}^{-1}b_2), \quad x_2 = b_2 - A_{21}x_1.$$
But it could happen that the Schur complement is not invertible. If $b_1 - A_{12}A_{22}^{-1}b_2$ is in the range of $[A/A_{22}]$ you have infinitely many solutions, since if $x_1^*$ is a solution, and $w$ is a non-zero vector from the kernel of $[A/A_{22}]$, then also $x_1^* + \alpha w$ is a solution.
Having one such solution $x_1^*$ you can compute $x_2^* = b_2 - A_{21}x_1^*$. In the paper they set $x_1^* = \frac{p_1+p_2}{2}$. I haven't proven that this is necessarily a solution for their matrix though.
The other case is that $b_1 - A_{12}A_{22}^{-1}b_2$ is
not in the range of $[A/A_{22}]$. Then there are no solutions. But you can try for the closest thing to a solution in an $L_2$ sense. Then $$x_1^+ = [A/A_{22}]^+(b_1 - A_{12}A_{22}^{-1}b_2),$$
where $[A/A_{22}]$ is the Moore-Penrose inverse. You may again set $x_2 = b_2 - A_{21}x_1^+$.
