This paper presents a new computational method for anisotropic tetrahedral meshing that (1) can control shapes of the elements by an arbitrary anisotropy function, and (2) can avoid ill-shaped elements induced from poorly distributed node locations. Our method creates a tetrahedral mesh in two steps. First our method obtains node locations through a physically based particle simulation, which we call 'bubble packing.' Ellipsoidal bubbles are closely packed on the boundary and inside a geometric domain, and nodes are placed at the centers of the bubbles. Our method then connects the nodes to create a tet mesh by the advancing front method. Experimental results show that our method can create a high quality anisotropic tetrahedral mesh that conforms well to the input anisotropy.