Using $20.3~{\rm fb}^{-1}$ of $e^{+}{e}^{-}$ collision data collected with the BESIII detector at $\sqrt{s}=3.773~{\rm GeV}$, we perform the first amplitude analysis of the decay $D^{+}\to {\pi }^{+}\eta \eta$. The intermediate process $D^{+}\to a_0(980{)}^{+}\eta$, $a_0(980{)}^{+}\to {\pi }^{+}\eta$, is the only component retained in the nominal amplitude model, and its branching fraction is measured to be $(3.67\pm0.12_{\rm stat}\pm0.06_{\rm syst})\times10^{-3}$. The ${\pi }^{+}\eta$ mass spectrum associated with $a_0(980{)}^{+}\eta$ production exhibits a line shape that differs substantially from those observed in $D_{(s)}\to a_0(980)\pi$ and $D^0\to a_0(980{)}^{-}{e}^{+}{\nu }_e$ decays. We examine several conventional descriptions of the $a_0(980)$ amplitude, including Flatt'e, dispersively modified Flatt'e, $T$-matrix, and $K$-matrix parameterizations. With reference $a_0(980)$ parameters, neither these models nor their extensions including additional small resonant or non-resonant amplitudes reproduce the observed line shape satisfactorily. When the $a_0(980)$ parameters are allowed to float, satisfactory fits can be obtained, but the pole mass is driven well above the $K\bar{K}$ threshold, inconsistent with the near-threshold character of the $a_0(980)$. The results reveal a tension between fit quality and the physical pole position in conventional direct-production amplitude models.