Concentrated multiarm star polymer solutions are used in this work as prototype, well defined, systems that clearly exhibit viscoplasticity and thixotropy. A thermodynamically consistent approach is used here to develop a continuum rheological model for these systems that is novel for thixotropic systems and has considerable advantages over previously postulated models: The basic rheological model is constructed from the beginning in a materially objective and flow type-independent form, in contrast to previous continuum models for thixotropy that typically are suitable only for shear flows. It is also consistent with non-equilibrium thermodynamics. The model is based on a Johnson-Segalman viscoelastic constitutive equation with a variable non-affine parameter. The non-affine parameter is considered representative of the structure and the aggregations between the particles in the system: as the deformation rate increases, those aggregations are assumed to be destroyed more and more, thus leading to discontinuities in the uniform viscoelastic deformation of the constituent multiarm polymer particles. This is described phenomenologically in the proposed model by a consistent increase in the non-affine parameter value from zero, which was taken to represent the virgin structure (where the deformation was assumed to be perfectly affine), to a maximum value at high rates, representing a concentrated suspension of viscoelastic particles. Simultaneously to those changes in the non-affine parameter, which are described by an empirical kinetic equation, the characteristic relaxation time of the system is also assumed to change dramatically: from infinity, at static equilibrium, to a minimum value characterizing the residual viscoelasticity of the fully broken-down medium. This behavior is modeled by assuming the relaxation time to be inversely proportional to the non-affine parameter. In addition to this viscoelastic contribution to the stress, a purely viscous contribution is considered modeled by a standard power law expression of the rate of strain. The model has been applied to transient and steady simple shear flows. The results show good semi-quantitative agreement with the experiments performed with the model multiarm star polymers.