A methodology for computing expansion basis functions using discrete harmonic modes is presented. The discrete harmonic modes are determined grain-by-grain for virtual polycrystals for which finite element meshes are available. The expansion weights associated with representing field variables over grain domains are determined by exploiting the orthogonality of the harmonic modes. The methodology is demonstrated with the representation of the axial stress distributions during tensile loading of a polycrystalline sample. An open source code, MechMonics, is available to researchers wishing to use the methodology to analyze data.