Setup
In this answer we'll use a concrete model as an illustration. It doesn't limit the generality of this discussion, but (hopefully) helps understanding.
In particular, we'll have a $z$-homogeneous (thus effectively 2D) non-magnetic $(\mu=1)$ crystal lit by a TE-polarized wave. The EM wave will then be described by a single equation for the $E_z$ component of the electric field, which we'll denote as simply $E$. We'll use the units where vacuum speed of light $c=1.$ The equation is then
$$\nabla^2 E=\varepsilon(\vec r)\partial^2_t E.\tag1$$
For a monochromatic wave with frequency $\omega$ we can write a time-independent equation:
$$\nabla^2 E=-\varepsilon(\vec r)\omega^2 E.\tag2$$
The concrete relative permittivity function $\varepsilon(\vec r)$ we'll take is
$$\varepsilon(\vec r)=1+(\varepsilon_{\text{max}}-1)\left(\sin\left(\frac{\pi x}a\right)\sin\left(\frac{\pi y}a\right)\right)^4.\tag3$$
We set the lattice period $a=1$ and peak permittivity $\varepsilon_{\text{max}}=3.5$.
The incident wave will come from the left at $\varphi=10°$ from the normal.
Symmetry and conserved quantities
Let's model our vacuum-crystal system the following way. Take the relative permittivity $\varepsilon(\vec r)$ from the crystal, aligning a crystal symmetry plane with $yOz$ coordinate plane. This crystal has a particular band structure with dispersion relations being $\omega_n(\vec K),$ where $\vec K$ is the quasi-wavevector, and $n$ is the number of band.
Now replace all the values of $\varepsilon(\vec r)$ for all $x<0$ with $1$, the vacuum value. For our example crystal $(3)$ the resulting $\varepsilon(x,y)$ will look as follows.
The replacement of the LHS of the crystal with vacuum breaks the original translation symmetry of the crystal. But partial symmetry—along $y$ and $z$ directions—still remains: any translation by a lattice vector lying in $yOz$ plane will yield the same $\varepsilon(\vec r)$. This means that the corresponding components of $\vec K$, i.e. $K_y$ and $K_z$, are still conserved on transition from the vacuum to the crystal.
For our example crystal, the unit strip (of the translation symmetry in the $y$ direction) for the resulting system will look as follows.
Reflected and transmitted modes
The remaining symmetry that we have lets us represent our full wavefunction as
$$E(\vec r)=u_{K_y,K_z}(\vec r)\exp(i(K_y y + K_z z)),\tag4$$
where $u_{K_y,K_z}$ is periodic in $y$ and $z$ per the Bloch theorem.
This means that to find the transmitted wave we only have to determine the $K_x$ component of the quasi-wavevector. To do this, we need to compute the band structure of the crystal—in particular, its cross section at the given values of $K_y$ and $K_z$.
Let's consider a wave incident on the crystal surface. It is described as
$$E_{\text{inc}}(\vec r)=E_0\exp(i\vec k\vec r),\tag5$$
$$\vec k=\omega\begin{pmatrix}\cos(\varphi)\\ \sin(\varphi)\\ 0\end{pmatrix}.\tag6$$
Taking $\omega=4.45$ and our values described above, we have
$$\vec k=\begin{pmatrix}4.3823945\\ 0.7727344\\ 0\end{pmatrix}.\tag7$$
The quasi-wavevector of the incident wave will then be
$$\vec K=\begin{pmatrix}-1.9007908\\ 0.7727344\\ 0\end{pmatrix}.\tag7$$
The band structure of the empty lattice (i.e. vacuum) with our example values at the slice of fixed $K_y=0.7727344$ will then look as follows.
From the intersection of $\omega$ shown by the gray line with the dispersion relation curves we can find the quasi-wavevector of the reflected wave. Only half of all the intersections correspond to the reflection, another half represent waves that propagate towards the crystal (one of which is the incident wave, others have zero amplitude).
Note how there are generally multiple possible reflected waves (e.g. at $\omega=8$). This is the result of diffraction on the periodic interface. If we do some calculations, we'll find that angles of this diffraction exactly follow the diffraction grating formula for different diffraction orders $m$:
$$\varphi_m^{\text{reflect}}=\pi-\arcsin\left(\sin\varphi-\frac{m\lambda}a\right).\tag8$$
Now, we are actually interested in the transmitted wave, rather than the reflected one. For this we should examine the band structure of the crystal at the same slice of $K_y.$
As we can see, generally there may be multiple transmitted waves. In our case there are only two intersections: one corresponds to the wave propagating to the right, and another to the left.
Direction of propagation of a crystal mode
Although choosing a particular crystal mode will let us find the phase velocity (though it may appear not so easy to identify for very high bands), this will tell nothing about the direction of a beam of light refracted by the crystal. This spatially localized behavior is described by group velocity instead.
Group velocity in a crystal is defined as usual:
$$\vec v_{\text{g}}=\nabla_{\vec K}\omega(\vec K).\tag9$$
Group velocity is also what we should use to choose from a pair of opposite $K_x$ to identify the transmitted wave: the transmitted mode is the one whose group velocity is directed into the crystal (i.e. in our example its $x$-component should be positive).
Snell's law and refractive index
Snell's law follows form conservation of $K_y$ and $K_z$ on transition through the interface. Considering the incident wave's quasi-wavevector in the incidence plane, we can write it as
$$\vec K_{\text{inc}}=K_{\text{inc}}\begin{pmatrix}\cos\varphi_{\text{inc}}\\ \sin\varphi_{\text{inc}}\end{pmatrix}.\tag{10}$$
Similarly for the transmitted wave we have
$$\vec K_{\text{trans}}=K_{\text{trans}}\begin{pmatrix}\cos\varphi_{\text{trans}}\\ \sin\varphi_{\text{trans}}\end{pmatrix}.\tag{11}$$
Since the $y$-components of $(10)$ and $(11)$ are equal, we get:
$$K_{\text{inc}}\sin\varphi_{\text{inc}}=K_{\text{trans}}\sin\varphi_{\text{trans}}.\tag{12}$$
If the material is homogeneous, then its Brillouin zone is infinite, and $\vec k=\vec K,$ and we recover the usual Snell's law by dividing both sides of $(12)$ by $k_{\text{inc}}.$ Otherwise, we have the generalized version of the Snell's law, where the refractive index is defined in terms of quasi-wavenumbers as
$$n=\frac{K_{\text{trans}}}{K_{\text{inc}}}.\tag{13}$$
Don't forget though, that crystals are generally anisotropic, so refractive index doesn't always let one find the angle of refraction of a beam of light. It's only near the center of the Brillouin zone in the lowest band that we can definitely say that $\omega(\vec K)$ becomes isotropic, so we can use Snell's law for this.
In some cases refractive index may become negative (our example is one such case), which will lead to an unusual refraction behavior.
