This research developed an automatic two-dimensional finite element meshing system to resolve practical engineering problems in the fields of geology, hydrology, and water resources. This system first used the Delaunay triangulation method to create reasonable-density triangular mesh and then converted it into quadrilateral mesh by combining proper pairs of adjacent triangles. A series of combination patterns aiming at three cases were established. The effect of the number of boundary edges on the subsequent meshing procedures were studied and summarized. For the geometry with multiple domains an adjustment method is proposed to completely eliminate the residual triangles during quadrilateral meshing through adjusting the number of boundary edges in each loop to be even. A special boundary loop identification method is proposed for priority treatment. Corresponding treatment methods aimed at three different situations are established for common boundary loops. For a certain boundary loop with an odd number of boundary edges, the appropriate edge for new point insertion is determined by the position properties and relative density errors. Practical applications confirm that the method proposed in this paper could successfully implement the full conversion from the triangular mesh to the quadrilateral mesh.