---

Abstract The element free Galerkin (EFG) method, which is based on the moving least square (MLS) approximation, requires only nodal data and no element connectivity. These features make the method more flexible than the conventional FEM. Nevertheless, direct imposition of the essential boundary conditions in the EFG method is always difficult because the shape functions obtained from the MLS approximation do not have the Kronocker-delta property. A new method named "the complementary integral method" is proposed here to overcome this difficulty. The presented method is more consistent with the variational basis of the EFG method. Several numerical examples are used to illustrate the implementation and performance of the method. The numerical examples including the Poisson's equation and 2D static and dynamic elasticity problems show that the method converges fast with reasonably accurate result for both the unknown variables and its derivatives.