Bubble packing method is extended to the case of three-dimensional domain to lead to an automatic tetrahedral finite element mesh. A static-equilibrium configuration is obtained from dynamic simulation for a system of bubbles packing a given domain. For this simulation, point mass for each bubble and viscous damping are taken into account together with reaction forces among the neighboring bubbles. The bubble population control is implemented to ensure a necessary number of bubbles. Discussed are the advantages of this physically based approach using a bubble located on a node. Comparison is made to other mesh generation schemes including the modified advancing front method.