This paper presents a computational method for quadrilateral meshing of a thin, or narrow, two-dimensional domain for finite element analysis. The proposed method creates a well-shaped single-layered, multi-layered, or partially multi-layered quadrilateral mesh. Element sizes can be uniform or graded. A high quality, layered quadrilateral mesh is often required for finite element analysis of a narrow two-dimensional domain with a large deformation such as in the analysis of rubber deformation or sheet metal forming. Fully automated quadrilateral meshing is performed in two stages: (1) extraction of the skeleton of a given domain by discrete chordal axis transformation, and (2) discretization of the chordal axis into a set of line segments and conversion of each of the line segments to a single quadrilateral element or multiple layers of quadrilateral elements. In each step a physically-based computational method called bubble packing is applied to discretize a curve into a set of line segments of specified sizes. Experiments show that the accuracy of a large-deformation FEM analysis can be significantly improved by using a well-shaped quadrilateral mesh created by the proposed method.