Yes, you are right, the conditional distribution needs to be found analytically, but I think there are lots of examples where the full conditional distribution is easy to find, and has a far simpler form than the joint distribution.
The intuition for this is as follows, in most "realistic" joint distributions $P(X_1,\dots,X_n)$, most of the $X_i$'s are generally conditionally independent of most of the other random variables. That is to say, some of the variables have local interactions, say $X_i$ depends on $X_{i-1}$ and $X_{i+1}$, but doesn't interact with everything, hence the conditional distributions should simplify significantly as $Pr(X_i|X_1, \dots, X_i) = Pr(X_i|X_{i-1}, X_{i+1})$
Ok, what you need to do is compute the joint posterior up to a constant, i.e. $f(y_1,...,y_n|\mu,\tau)p(\mu,\tau)$. Then to compute the conditional posterior $\pi(\mu|\mathbf{y},\tau)$ you just treat the $\tau$ terms as fixed and known, so that some of them can be cancelled out. Then you do the same thing with $\mu$ to get $\pi(\tau|\mathbf{y},\mu)$.
So what I mean is, the joint posterior is given by:
\begin{equation}
\pi(\mu,\tau|\mathbf{y}) \propto \tau^{\frac{n}{2}+\alpha_0-1} \exp\{-\frac{\tau}{2} \sum_{i=1}^{n}(y_i-\mu)^2-\frac{\tau_0}{2}(\mu-\mu_0)^2-\beta_0\tau\}.
\end{equation}
Now, to get the conditional posterior $\pi(\tau|\mathbf{y},\mu)$ you just remove the $\mu$ terms that just multiply the expression. So we could re-write the joint posterior as:
\begin{equation}
\pi(\mu,\tau|\mathbf{y}) \propto \tau^{\frac{n}{2}+\alpha_0-1} \exp\{-\frac{\tau_0}{2}(\mu-\mu_0)^2\}\exp\{-\frac{\tau}{2} \sum_{i=1}^{n}(y_i-\mu)^2-\beta_0\tau\}.
\end{equation}
Now, since we are treating $\mu$ as fixed and known, for the conditonal posterior $\pi(\tau|\mu,\mathbf{y})$ we can remove the first exponential term from the equation and it will still be true (owing to the $\propto$ sign rather than the equals sign). It's important that you get this: the conditional posterior says given that we know $\mu$ what is $\tau$, so $\mu$ is known (and hence fixed). So the kernel for the conditional posterior for $\tau$ becomes (after some re-arranging):
\begin{equation}
\pi(\tau|\mu,\mathbf{y}) \propto \tau^{\frac{n}{2}+\alpha_0-1}\exp\{-\tau(\frac{1}{2} \sum_{i=1}^{n}(y_i-\mu)^2+\beta_0)\},
\end{equation}
which is the kernel for the Gamma distribution with the parameters you state.
Now, for $\pi(\mu|\tau,\mathbf{y})$ you do the same thing, but it's a bit trickier because you have to complete the square. After cancelling the relevant terms not involving $\mu$ you get:
\begin{equation}
\pi(\mu|\tau,\mathbf{y}) \propto \exp\{ -\frac{\tau}{2} \sum_{i=1}^{n} (y_i-\mu)^2 - \frac{\tau_0}{2}(\mu-\mu_0)^2\}.
\end{equation}
If you multiply out the squared brackets and take the relevant summations, then remove all terms not involving $\mu$ (as they are all just multiplying the kernel by a constant) this becomes:
\begin{equation}
\pi(\mu|\tau,\mathbf{y}) \propto \exp\{ -\frac{1}{2}[\mu^2(n\tau+\tau_0) - 2\mu(\tau n\bar{y} - \tau_0\mu_0)]\}.
\end{equation}
If you complete the square here you get the kernel for a Gaussian with the mean and variance you state.
Hope that helps.
Best Answer
Reading through the comments on the othe answers, I believe the correct answer to the question that was intended to be asked is "they don't", in general. As has been mentioned, they construct a DAG and look at the Markov blanket and then (roughly) do the following.
This isn't exactly what is being done; for example, JAGS will use some other tricks to construct block updates. But this should give an idea of what they are doing.