Assume we are on a smooth, complete Riemannian manifold $(M,g), dim(M) \geq 3$. What are the specific geometric/topological constraints for such a manifold to admit complete, totally geodesic hypersurfaces? Please, I admit that I'm a rookie so any simple illustrations and sources are more than welcome, especially for the case of lowest dimension $(dim(M) = 3)$. 



existence of totally geodesic hypersurfaces
You should be aware that, for $n\ge3$, the generic Riemannian metric $(M^n,g)$ has no totally geodesic hypersurfaces at all, even locally. Typically, when they do exist, it is for some geometric reason that makes them obvious. For example, if $(M^n,g)$ admits an isometry $\iota: M\to M$ of order $2$ (i.e., an involution) whose fixed point set is a hypersurface $H^{n-1}\subset M$, then $H$ will be totally geodesic (and it will be complete if $g$ is complete).


As for topological constraints, there really aren't any, because, for any closed, embedded hypersurface $H^{n-1}\subset M^n$, you can always construct a metric $g$ such that $H$ is totally geodesic in $(M,g)$; just take the metric to be a product metric on some tubular neighborhood of $H$.


As a practical matter, if someone hands you a Riemannian manifold $(M^n,g)$ that is described sufficiently explicitly that you can compute the various curvature tensors and functions to arbitrary orders and determine their vanishing loci, there is an algorithm to follow that will determine all of the totally geodesic hypersurfaces in $(M^n,g)$, and that is probably the best you can do in terms of giving explicit 'constraints' on $(M^n,g)$ for it to admit totally geodesic hypersurfaces (complete or not). The algorithm basically, goes like this: Let $\pi: \Sigma^{2n-1}\to M^n$ be the unit sphere bundle of $M$ with respect to the metric $g$. Then, using the Levi-Civita connection, construct an $(n{-}1)$-plane field $D\subset T\Sigma$ such that a curve $\nu:\mathbb{R}\to \Sigma$ of the form $\nu(t) = \bigl (x(t), e(t)\bigr)$, where $e(t)\in T\_{x(t)}M$ is a unit vector, is tangent to $D$ if and only if $x'(t)\cdot e(t)=0$ and $e(t)$ is parallel along $x(t)$, i.e., $\nabla\_{x'(t)}e(t) = 0$. Then a hypersurface $H\subset M$ is totally geodesic if and only if its 'unit normal graph' in $\Sigma$ is everywhere tangent to $D$. (It is easy to show that this plane field $D$ is Frobenius if and only if $(M^n,g)$ has constant sectional curvature, which explains why only those manifolds have totally geodesic hypersurfaces through every point perpendicular to every direction.) Now the problem is reduced to finding the $(n{-}1)$-manifolds in $\Sigma$ that are tangent to $D$, which is a standard problem.


It is probably worth summarizing what this process says about Riemannian $3$-manifolds $(M^3,g)$. When one goes through the above analysis, the first thing it tells one is that, if $H^2\subset M$ is a totally geodesic hypersurface with normal vector field $\nu:H\to \Sigma$, then, for all $x\in H^2$, the normal vector $\nu(x)\in T\_xM$ must be an eigenvector of the Ricci curvature $\mathrm{Ric}(g)$. In particular, if $\mathrm{Ric}(g)$ has distinct eigenvalues everywhere (or, what would be generic, on a dense open set $U\subset M$), then one can locally write
$$
g = {\omega\_1}^2+{\omega\_2}^2+{\omega\_3}^2
\quad\text{and}\quad
\mathrm{Ric}(g) = \lambda\_1\,{\omega\_1}^2+ \lambda\_2\,{\omega\_2}^2+ \lambda\_3\,{\omega\_3}^2
$$
for some functions $\lambda\_1 < \lambda\_2 < \lambda\_3$ and some coframing $\omega\_i$, 
and one of the $\omega\_i$ will have to vanish on $H$. However, if $\omega\_i$ vanishes on $H$, then $\mathrm{d}\omega\_i$ must also vanish on $H$. Defining the functions $f\_{ij}$ that satisfy $\omega\_i\wedge\mathrm{d}\omega\_j = f\_{ij}\,\omega\_1\wedge\omega\_2\wedge\omega\_3$, one sees that $f\_{ii}$ must also vanish on $H$. 


In the generic situation, each of the three equations $f\_{ii} = 0$ will define a surface $H\_i$ in $M$, and these $H\_i$ are the only possible totally geodesic surfaces in $M$. In fact, $H\_i$ will be totally geodesic if and only if $\omega\_i$ vanishes on $H\_i$ along with the three functions $f\_{jk}$, $f\_{kj}$ and $f\_{jj}{-}f\_{kk}$, where $(i,j,k)$ is a permutation of $(1,2,3)$. Thus, this provides the effective test in the case that the Ricci tensor has distinct eigenvalues. 


In the degenerate case of multiple eigenvalues (or, say, that $f\_{ii}=0$ does not define a surface in $M$), more differentiation is required to formulate an explicit test of this nature, but I'll leave that to you.


