  We propose a nonlocal operator method for solving partial differential
equations (PDEs). The nonlocal operator is derived from the Taylor series
expansion of the unknown field, and can be regarded as the integral form
"equivalent" to the differential form in the sense of nonlocal interaction. The
variation of a nonlocal operator is similar to the derivative of shape function
in meshless and finite element methods, thus circumvents difficulty in the
calculation of shape function and its derivatives. {The nonlocal operator
method is consistent with the variational principle and the weighted residual
method, based on which the residual and the tangent stiffness matrix can be
obtained with ease.} The nonlocal operator method is equipped with an hourglass
energy functional to satisfy the linear consistency of the field. Higher order
nonlocal operators and higher order hourglass energy functional are
generalized. The functional based on the nonlocal operator converts the
construction of residual and stiffness matrix into a series of matrix
multiplications on the nonlocal operators. The nonlocal strong forms of
different functionals can be obtained easily via support and dual-support, two
basic concepts introduced in the paper. Several numerical examples are
presented to validate the method.
