This paper presents a new polyomino tiling algorithm for three-dimensional surface geometric model. We first apply surface subdivision to the input polygonal mesh or subdivision surface mesh. Then we construct a random Hamiltonian path to connect all subdivided input mesh. This path is used to build initial polyomino tiling on the original input mesh. Finally we apply random polyomino exchanging to the polyomino tiling to get more uniform occurrence of each polyomino types. Our method is applicable to construct three-dimensional puzzle and we show the results of proposed algorithm on three-dimensional mesh data.