Coupled thermo-hydro-mechanical (THM) processes are ubiquitous in subsurface energy production and geological utilization and storage operations. Numerical simulation of strongly coupled THM processes is a non-trivial task, yet required to predict the performance of many applications in energy geomechanics. The majority of existing and open source THM numerical codes are not end user adaptable and do not include elastoplasticity coupled to mass and energy balance equations. This article presents an open source thermo-poroelastoplastic finite element numerical code with a fully-coupled monolithic solution strategy that is solved with Fenicsx computing platform. The formulation employs a mixed finite element scheme for pore pressure diffusivity, Petrov-Galerkin methods for energy transport, mean stress dependent yield surface, and non-associative plastic potential. The numerical solution is verified with small-scale conventional triaxial tests, including drained and undrained compression and extension. We present example simulations reaching the yield surface induced by coupled hydro-mechanical and thermal loads. In addition, we present two example large-scale applications related to geothermal energy and carbon geological storage. Results show that the numerical solution accurately predicts changes of temperature, pore pressure, and stress for a wide range of model geometries and boundary conditions, including the plastic response. The code is freely available to the general community for use and modification.