Abstract: By absorbing advantages of the finite element and meshless methods, a new numerical method, free element method, is proposed in the paper. In the discretization, the isoparametric elements as used in FEM are employed to represent the geometry and interpolate physical variables; and in the algorithm, the point collocation technique using elements is employed to generate the system of equations point by point. The main feature of the method is that only one independent element formed by freely selecting surrounding points is required for each collocation point, without need to consider the connective relationship between adjacent elements and the continuity of physical variables and their spatial derivatives at interfaces of the connected elements. Two types of free element methods, the strong-form method and weak form method, will be described in the paper. The former directly generates the system of equations from the governing equations and the Neumann boundary conditions, while the latter establishes the weak-form integral expression of the governing equations by the weighted residual technique over the free element first and then generates the system of equations through an integration process similar to that employed in the standard FEM. The method proposed in the paper is an element collocation method. To achieve highly accurate spatial derivatives for internal collocation points of the computational domain, isoparametric elements with at least one internal node are required. For this purpose, apart from the arbitrary order quadrilateral Lagrange elements, a new seven-node triangle element is constructed in the paper, which can be used to model problems with complex geometries.