Interpretation of a singular metric
The nature of singularities in GR is a delicate issue. A good review of the difficulties presented to define a singularity are in Geroch's paper What is a singularity in GR?
The problem of attaching a boundary in general to a spacetime is that there is not natural way to do it. for example, in the FRW metric the manifold at $t=0$ can be described by two different coordinate systems as: $$\{t,r\cos\theta,r\sin\theta \cos\phi,r\sin\theta \sin\phi\}$$ or $$\{t,a(t)r\cos\theta,a(t)r\sin\theta \cos\phi,a(t)r\sin\theta \sin\phi\}$$ In the first case we have a three dimensional surface, in the latter a point.
It might be tempting to define a singularity following other physical theories as the points where the metric tensor is undefined or below $C^{2}$. However, this is troublesome because in the gravitational case the field defines also the spacetime background. This represents a problem because the size, location and shape of singularities can't be straightforward characterize by any physical measurement.
The theorems of Hawking and Penrose, commonly used to show that singularities in GR are generic under certain circumstances have the conclusion that spacetime must be geodesically incomplete (Some light-paths or particle-paths cannot be extended beyond a certain proper-time or affine-parameter).
As mentioned above the peculiar characteristic of GR of identifying the field and the background makes the task of assigning a location, shape or size to the singularities very delicate. If one thinks in a singularity of the gravitational potential in classical terms the statement that the field diverges at a certain location is unambiguous. As an example, take the gravitational potential of a spherical mass $$V(t,r,\theta,\phi)=\frac{GM}{r}$$ with a singularity at the point $r=0$ for any time $t$ in $\mathbb{R}$. The location of the singularity is well defined because the coordinates have an intrinsic character which is independent of $V$ and are defined with respect the static spacetime background.
However, this prescription doesn't work in GR. Consider the spacetime with metric $$ds^{2}=-\frac{1}{t^{2}}dt^{2}+dx^{2}+dy^{2}+dz^{2}.$$ defined on $\{(t,x,y,z)\in \mathbb{R}\backslash \{0\}\times \mathbb{R}^{3}\}$. If we say that there is a singularity at the point $t=0$ we might be speaking to soon for two reasons. The first is that $t=0$ is not covered by our coordinate chart. It makes no sense to talk about $t=0$ as a point in our manifold using these coordinates. The second thing is that the lack of an intrinsic meaning of the coordinates in GR must be taken seriously. By making the coordinate transformation $\tau=\log(t)$ we obtain the metric $$ds^{2}=d\tau^{2}+dx^{2}+dy^{2}+dz^{2},$$ on $\mathbb{R}^{4}$ and remain isometric to the previous spacetime defined in $\{(t,x,y,z)\in \mathbb{R}\backslash \{0\}\times \mathbb{R}^{3}\}$. What we have done is find an extension of the metric to $\mathbb{R}^{4}$. The singularity was just a coordinate singularity, similar to the event horizon singularity in Schwarzschild coordinates. The extended spacetime is of course Minkowski spacetime which is non-singular.
Another approach is to define a singularity in terms of invariant quantities such as scalar polynomials of the curvature.
This are scalars formed by the Riemann tensor. If this quantities diverge it matches our physical idea that and object approaching regions of higher and higher values must suffer stronger and stronger deformations. Also, in many relevant cosmological models like FRW and Black Holes metrics one can show that this indeed happen. But as mentioned the domain of the gravitational field defines the location of events so a point where the curvature blow up might not be even in the domain. Therefore, we must formalise the following: statement "The scalar diverges as we approach a point that has been cut out of the manifold.". If we were in a Riemann manifold then the metric define a distance function $$d(x,y):(x,y)\in\cal{M}\times\cal{M}\rightarrow \inf\left\{\int\rVert\dot{\gamma}\rVert \right\}\in\mathbb{R}$$
where the infimum is taken over all piecewise $C^{1}$ curves $\gamma$ from $x$ to $y$. Moreover, the distance function allows us to define a topology. A basis of that topology is given by the set $\{B(x,r):y\in{\cal{M}}| d(x,y)\le r \forall x\in \cal{M}\}$. The topology naturally induce a notion of convergence. We say the sequence $\{x_{n}\}$ converges to $y$ if for $\epsilon> 0$ there is an $N\in \mathbb{N}$ such that for any $n\ge N$ $d(x_{n},y)\le \epsilon$. A sequence that satisfies this conditions is called a Cauchy sequence. If every Cauchy sequence converges we say that $\cal{M}$ is metrically complete
Notice that now we can describe points that are not in the manifold as a point of convergence of a sequence of points that are. Then the formal statement can be stated as: "The sequence $\{R(x_{n})\}$ diverges as the sequence $\{x_{n}\}$ converges to $y$" where $R(x_{n})$ is some scalar evaluated at $x_{n}$ in $\cal{M}$ and $y$ is some point not necessarily in $\cal{M}$.
In the Riemannian case if every Cauchy sequence converges in $\cal{M}$ then every geodesic can be extend indefinitely. That means we can take as the domain of every geodesic to be $\mathbb{R}$. In this case we say that $\cal{M}$ is geodesically complete. In fact also the converse is true, that is if $\cal{M}$ is geodesically complete then $\cal{M}$ is metrically complete.
So far, all the discussion has been for Riemann metrics, but as soon as we move to Lorentzian metrics the previous discussion can't be used as stated. The reason is that Lorentzian metrics doesn't define a distance function. They do not satisfy the triangle inequality. So we only have left the notion of geodesic completeness.
The three kinds of vectors available in any Lorentzian metric define three nonequivalent notions of geodesic completeness depending on the character of the tangent vector of the curve: spacelike completeness, null completeness and timelike completeness. Unfortunately, they are not equivalent it is possible to construct spacetimes with the following characteristics:
- timelike complete, spacelike and null incomplete
- spacelike complete, timelike and null incomplete
- spacelike complete, timelike and null incomplete
- null complete, timelike and spacelike incomplete
- timelike and null complete, spacelike incomplete
- spacelike and null complete, timelike incomplete
- timelike and spacelike complete, null incomplete
Moreover, in the Riemannian case if $\cal{M}$ is geodesically complete it implies that every curve is complete, that means every curve can be arbitrarily extended . Again, in the Lorentzian case that is not the case, Geroch construct an example of a geodesically null, timelike and spacelike complete spacetime with a inextendible timelike curve of finite length. A free falling particle following this trajectory will accelerate but in a finite amount of time its spacetime location would stop being represented as a point in the manifold.
Schmidt provided an elegant way to generalise the idea of affine length to all curves, geodesic and no geodesics. Moreover, the construction in case of incomplete curves allows to attach a topological boundary $\partial\cal{M}$ called the b-boundary to the spacetime $\cal{M}$.
The procedure consist in building a Riemannian metric in the frame bundle $\cal{LM}$. We will use the solder form $\theta$ and the connection form $\omega$ associated to the Levi-Civita connection $\nabla$ on $\cal{M}$. Explicitly,
\begin{equation} G_{ab}(X_{a},Y_{a})= \theta(X_{a}) \cdot \theta(Y_{a})+\omega(X_{a})\bullet \omega(Y_{a}) \end{equation}
where $X_{a},Y_{a}\in T_{p}\cal{LM}$ and $\cdot,\bullet$ are the standard inner product in $\mathbb{R}^{n}$ and $\mathfrak{g}\cong\mathbb{R}^{n^{2}}$.
Let $\gamma$ be a $C^{1}$ curve through $p$ in $\cal{M}$ and a basis $\{E_{a}\}$. Now choose a point $P$ in $\cal{LM}$ such that $P$ satisfies $\pi(P)=p$ and the basis of $T_{p}$ is given by $\{E_{a}\}$. Using the covariant derivative induced by the metric we can parallel propagate $\{E_{a}\}$ in the direction of $\dot{\gamma}$. This procedure defines a curve $\Upsilon $ in $\cal{LM}$. This curve is called the lift of the curve $\gamma$. The length of $\Upsilon$ with respect to Schmidt metric, $$l=\int_{\tau}\|\dot{\Upsilon} \|_{G} dt$$ is is called a generalised affine parameter. If $\gamma$ is a geodesic $l$ is an affine parameter. If every curve in a spacetime $\cal{M}$ that has finite affine generalised length has endpoints we call the spacetime b-complete. If it is not $b$- complete we call the spacetime b-incomplete. A classification of singularities in terms of the b-boundary (See Chapter 8, The large Scale Structure of spacetime) was done by Ellis and Schmidt here.
In the case of the FRW the b-boundary $\partial\cal{M}$ was computed in this paper The result is that the boundary is a point. However, the resulting topology in $\partial\cal{M}\cup\cal{M}$ is non-Hausdorff. This means the singularity is in some sense arbitrary close to any event in spacetime. This was regarded as unphysical and attempts to improve the b-boundary construction were made without any attempt having a particular acceptance. Also the high dimensionality of the bundles involved make the b-boundary a difficult working tool.
Another types of boundaries can be attached. For example:
conformal boundaries used in Penrose diagrams and in the AdS/Cft correspondance. In this case the conformal boundary as seen here at $t=0$ is a three dimensional manifold.
Causal boundaries. This constructions depends only on the causal structure, so it doesn't distinguish between boundary points at a finite distance or at infinity. (See chapter 6, The large scale structure of spacetime)
Abstract boundary.
I am unaware if in the two last cases explicit calculations have been done to the case of the FRW metric.
Just for clarification: the manifolds used in general relativity are locally (in the sense of diffeomorphisms) $\mathbb R^{3+1}$. They are not $ \mathbb R^{3+1}$ in general, which is Minkowski space with zero curvature.
In this sense, the points at $t=0$ do not belong to the manifold, as there is no neighborhood which is diffeomorphic to $\mathbb R^{3+1}$ as @Ali Moh already stated. This means that with the FLRW metric alone, one cannot make predictions for the "big bang" (although one can make predictions around $t=0+\epsilon$ for every $\epsilon > 0$).
I'm assuming you are not talking about the modern view of the history of the universe where the big bang (reheating) was preceded by inflation, and we don't know what happened before inflation or whether there is even a beginning of time.
So if you want to just look at the FRW metric then I think you're comment that the singularity does not belong to the manifold is correct because the limit $t \rightarrow 0$ is not well defined. For example looking at the volume $$\lim_{t\rightarrow 0^+} \int d^4 \sqrt{-g} = \infty$$ whereas $$\int d^4 \sqrt{-g}\,\,\Big|_{t=0} = 0$$
So you cannot smoothly add the singularity point to the manifold.