This is a typical question in finite element approximation for incompressible Stokesian flow. Also since you mentioned the paper of Crouzeix and Thomée, both of whom are very famous in finite element community, a finite element counter-example here seems appropriate.
The question is:
Is the projection bounded? If we project $L^2$-regular divergence-free vector fields into a finite dimensional subspace of $H^1$-regular divergence-free vector fields.
The answer is No.
Spaces: Consider $\Omega\in \mathbb{R}^2$ being a simply-connected bounded smooth domain
$$\newcommand{\b}{\boldsymbol}
H = \{\b{v}\in L^2(\Omega): \mathrm{div}\,\b{v} = 0 \},
$$
and
$$
V = \{\b{v}\in H^1(\Omega): \mathrm{div}\,\b{v} = 0 \}.
$$
$V$ is then dense in $H$ (for a proof please refer to my answer here: Relation of the kernels of one bounded operator and its extension , basically I modified Tartar's argument in his book for $C^{\infty}$-vector fields, a much simpler argument would be simply using $C^{\infty}(\Omega)\subset H^1(\Omega)\subset L^2(\Omega)$ and density of $C^{\infty}$-vector fields).
Norms and continuous embedding: Norms for $H$ is the $H(\mathrm{div})$-norm :
$$
\|\cdot \|_{H(\mathrm{div})} := \left(\|\cdot \|_{L^2}^2 + \|\mathrm{div }\cdot \|_{L^2}^2\right)^{1/2}.
$$
Norm for $V$ is the standard $H^1$-norm:
$$
\|\cdot \|_{H^1} := \left(\|\cdot \|_{L^2}^2 + \|\nabla (\cdot) \|_{L^2}^2\right)^{1/2}.
$$
Clearly for every $\b{v}\in V$:
$$\left(\|\b{v} \|_{L^2}^2 + \|\mathrm{div }\,\b{v} \|_{L^2}^2\right)^{1/2} = \|\b{v}\|_{H}\leq \|\b{v}\|_{V}=\left(\|\b{v}\|_{L^2}^2 + \|\nabla \b{v}\|_{L^2}^2\right)^{1/2}$$ because the divergence's $L^2$-norm is bounded by the full gradient's $L^2$-norm, not just the incompressible ones (incompressible means divergence free), hence continuously embeddedness is true as well.
Going into finite dimension: Let $V_m \subset V$ be a finite element space, of which each element is a piecewise continuous polynomial vector field. Please see Thomasset's book in finite element for Navier-Stokes equation for a construction on a Delaunay triangulation, or Brenner and Scott's finite element book exercise 11.x.24, both of which have the formula for the polynomial vector fields.
Now consider the projection
$$P_m: H\to V_m\cap\{\b{v} = 0\text{ on }\partial \Omega\}\subset V$$ where the orthogonality is in $H$'s inner product: for any $\b{w}_i\in V_m$
$$
(P_m \b{v} - \b{v},\b{w}_i)_{L^2} + (\mathrm{div}(P_m \b{v} - \b{v}),\mathrm{div}\,\b{w}_i)_{L^2} = 0.
$$
Notice this is the same as orthogonality in $L^2$ because of divergence free:
$$
(P_m \b{v} - \b{v},\b{w}_i)_{L^2} =0.
$$
This only implies the projection is stable when measured under $L^2$-norm. Now rewrite $\b{v}$'s $H^1$-norm if it vanishes on the boundary:
$$
\|\b{w}\|^2_{H^1} = \|\b{w}\|^2_{L^2} + \|\mathrm{div}\,\b{w}\|^2_{L^2} + \|\mathrm{curl}\,\b{w}\|^2_{L^2} .
$$
Also notice $\b{v}$ is divergence free, therefore the stability we want to get is bounded by
$$
\frac{\|P_m \b{v}\|_{V}}{\|\b{v}\|_{H}} = \frac{\|P_m \b{v}\|_{H^1}}{\|\b{v}\|_{H(\mathrm{div})}} = \frac{\left(\|P_m \b{v}\|_{L^2}^2 + \|\mathrm{curl}(P_m \b{v})\|_{L^2}^2\right)^{1/2}}{\|\b{v}\|_{L^2 }}.\tag{1}
$$
Notice that the curl of the projection $\b{v}$ can be unbounded for $m\to \infty$, for the curl's $L^2$-norm can be written as the sum of the norm on each element $K$ in this triangulation
$$
\|\mathrm{curl}(P_m \b{v})\|_{L^2}^2 = \sum_{K\in \mathcal{T}}\|\mathrm{curl}(P_m \b{v})\|_{L^2(K)}^2
$$
Now use Poincaré inequality for divergence free vector fields with vanished boundary:
$$
\|\b{p} - \bar{\b{p}}|_M\|_{L^2(D)} \leq C(D) \|\mathrm{curl}\,\b{p}\|_{L^2(D)},
$$
where $\bar{\b{p}}_D = \frac{1}{|D|}\int_D \b{p} \,d\b{x}$ is the averaging vector field on $D$. and the Poincaré constant is the same order as the $n$-th root of measure of $D$ where $n$ is the dimension of the $D$, in our case the dimension is $n=2$. Therefore
$$
\sum_{K\in \mathcal{T}}\|\mathrm{curl}(P_m \b{v})\|_{L^2(K)}^2 = \sum_{K\in \mathcal{T}}\|\mathrm{curl}(P_m \b{v} - \overline{P_m \b{v}}_K)\|_{L^2(K)}^2 \\
\geq \sum_{K\in \mathcal{T}} O(|K|^{-1}) \| P_m \b{v} - \overline{P_m \b{v}}_K \|_{L^2(K)}^2
\\
\geq \min_{K\in \mathcal{T}} |K|^{-1} \sum_{K\in \mathcal{T}} \| P_m \b{v} - \overline{P_m \b{v}}_K \|_{L^2(K)}^2 \geq O(M)\| P_m \b{v} - \overline{P_m \b{v}} \|_{L^2 }^2,\tag{2}
$$
where $M = \dim V_m$, which is the number of element in this finite dimensional approximation space $V_m\subset V$. $\overline{P_m \b{v}}$ restricted on $K$ is the local averaging vector field.
To see why (1) may be unbounded, we can choose a highly oscillated sequence $\{\b{v}_m\}\subset H^1$ so that its oscillation can not be resolved by the density of this triangulation at level $m$, so that for every $V_m$:
$$
\| P_m \b{v}_m - \overline{P_m \b{v}_m} \|_{L^2 }\sim \| P_m \b{v}_m\|_{L^2 }.
$$
Think divergence free vector field sequence
$$
\b{v}_m = \big(\sin(Mx)\cos(My), -\cos(Mx)\sin(My)\big)\in H^1([0,\pi]^2),
$$
where $M$ is the number of the simplices (triangles in this case) in level $m$ triangulation $\mathcal{T}_m$.
Inequality in (2) is true for the dimension of $V_m$ depends on the number of the simplices in this triangulation of $\Omega$, so that if we say $|\Omega| = O(1)$, then $O(|K|^{-1})$ is the same order as $M$. If the triangulation goes finer and finer, $|K|\to 0$ for every $K\in \mathcal{T}$, $M\to \infty$, and (1) becomes unbounded for this sequence.
Backstory if you are interested in finite element for incompressible flow: using $H^1$-divergence free elements to approximate $H(\mathrm{div})$-regular vector fields is not a good idea based on this.
