The article offers a general formulation for the multipolar field expansion of an inclusion. The multipolar moments of arbitrary order for the ellipsoidal inclusion are derived, showing known features of their elastic fields in both the far and near field. Using integer point enumeration theory, the same formulation is extended to all polytopal inclusions, which serve to model faceted inclusions found in, e.g. igneous and metamorphic rocks, particle reinforced composite materials, or as intermetallic precipitates in alloys. A general algorithm for the calculation of the multipolar moments of polytopes of any number of vertices is derived. A number of examples for tetrahedral, cuboidal and dodecahedral inclusions are given. It is shown that the parity of the axial moment function of a polytopal inclusion mirrors the symmetry of the inclusion. This represents the main properties of their elastic fields, including the long- and near-field decay. If the inclusion is symmetric with respect to a given axes, the fields will decay with
1
/
r
2
in the far field and
1
/
r
4
in the near field. If symmetry along that direction is lost, as is expected in most real inclusions, the decay rate progresses with
1
/
r
2
but the near field with
1
/
r
3
.