def METHOD_NAME(self, n):
    """Uniform boundary points on the spatio-temporal domain.
    Geometry surface area ~ bbox.
    Time surface area ~ diam.
    """
    if self.geometry.dim == 1:
        nx = 2
    else:
        s = 2 * sum(
            map(
                lambda l: l[0] * l[1],
                itertools.combinations(
                    self.geometry.bbox[1] - self.geometry.bbox[0], 2
                ),
            )
        )
        nx = int((n * s / self.timedomain.diam) ** 0.5)
    nt = int(np.ceil(n / nx))
    x = self.geometry.METHOD_NAME(nx)
    nx = len(x)
    t = np.linspace(
        self.timedomain.t1,
        self.timedomain.t0,
        num=nt,
        endpoint=False,
        dtype=config.real(np),
    )
    xt = []
    for ti in t:
        xt.append(np.hstack((x, np.full([nx, 1], ti))))
    xt = np.vstack(xt)
    if n != len(xt):
        print(
            "Warning: {} points required, but {} points sampled.".format(n, len(xt))
        )
    return xt