Engineering problems frequently require solution of governing equations with spatially-varying discontinuous coefficients. Even for linear elliptic problems, mapping large ensembles of coefficient fields to solutions can become a major computational bottleneck using traditional numerical solvers. Furthermore, machine learning methods such as neural operators struggle to fit these maps due to sharp transitions and high contrast in the coefficient fields and a scarcity of informative training data. In this work, we focus on a canonical problem in computational mechanics: prediction of local elastic deformation fields over heterogeneous material structures subjected to periodic boundary conditions. We construct a hybrid approximation for the coefficient-to-solution map using a Thermodynamically-informed Iterative Neural Operator (TherINO). Rather than using coefficient fields as direct inputs and iterating over a learned latent space, we employ thermodynamic encodings -- drawn from the constitutive equations -- and iterate over the solution space itself. Through an extensive series of case studies, we elucidate the advantages of these design choices in terms of efficiency, accuracy, and flexibility. We also analyze the model's stability and extrapolation properties on out-of-distribution coefficient fields and demonstrate an improved speed-accuracy tradeoff for predicting elastic quantities of interest.