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^+$.