Abstract: The governing equations of strain gradient elasticity are fourth-order partial differential equations. Galerkin discretization would require C¹ continuous shape functions. C¹ natural neighbor interpolant can be realized when embedding non-Sibsonian interpolant in the Bernstein-Bézier surface representation of a cubic simplex. Essential boundary conditions are imposed directly in a Galerkin scheme for the strain gradient elasticity because the C¹ interpolant has the interpolating properties for nodal function and nodal gradient values. Boundary layer analysis and infinite plate with a central circular hole subjected to the biaxial tension are analyzed. The numerical solutions agree well with the analytical solutions, which shows that C1 natural neighbor Galerkin method can be used to analyze the strain gradient elasticity problems.