Hexagonal grid methods are found useful in many research works, including numerical modeling in spherical coordinates, in atmospheric and ocean models, and simulation of electrical wave phenomena in cardiac tissues. Almost all of these used standard Laplacian and mostly on one configuration of regular hexagons. In this work, discrete symmetric boundary condition and energy product for anisotropic Laplacian are investigated firstly on general net of regular hexagons, and then generalized to its most extent in two-or three-dimensional cell-center finite difference applications up to the usage of symmetric stencil in central differences. For analysis of Laplacian related applications, this provides with an approach in addition to the M-matrix theory, series method, functional interpolations and Fourier vectors.