What's the most general form of Maxwell's stress tensor for EM fields in matter?

In Griffith's Introduction to Electrodynamics we find a derivation of the relation $$\mathbf f + \frac{\partial \mathbf g}{\partial t} = \nabla \cdot \overleftrightarrow{\mathbf T} ,$$ where $$\mathbf f$$ is the Lorentz force per unit volume, $$\mathbf g$$ is the electromagnetic linear momentum per unit volume, and $$\overleftrightarrow{\mathbf T}$$ is Maxwell's stress tensor. However, this derivation employs Maxwell's equations in a vacuum, or at least in media that exhibit no polarization or magnetization.

My aim is to derive a similar relation for Maxwell's tensor, but starting from the more general set of equations that include the fields $$\mathbf D$$, $$\mathbf E$$, $$\mathbf B$$, and $$\mathbf H$$, with the general relations $$\mathbf D = \varepsilon_0 \mathbf E + \mathbf P$$ and $$\mathbf B = \mu_0(\mathbf H + \mathbf M)$$. Starting again from the local expression of the Lorentz force on free charges and current densities $$\mathbf f = \rho_0 \mathbf E + \mathbf J_0 \times\mathbf B,$$ I've managed to show that $$\mathbf f = -\frac{\partial \mathbf g}{\partial t} + (\nabla \cdot \mathbf D)\mathbf E + (\nabla \cdot \mathbf B) \mathbf H - \mathbf D \times (\nabla \times \mathbf E) - \mathbf B \times (\nabla \times \mathbf H),$$ but I cannot seem to be able to continue. Intuitively, I think I should get something of the form $$T_{ij} = E_i D_j + H_i B_j - \frac 1 2 \delta_{ij} (\mathbf E \cdot \mathbf D + \mathbf H \cdot \mathbf B),$$ which does revert to Griffith's $$T_{ij} = \varepsilon_0 E_i E_j + \frac 1 {\mu_0} B_i B_j - \frac 1 2 \delta_{ij} \left( \varepsilon_0 E^2 + \frac 1 {\mu_0} B^2\right)$$ as soon as $$\mathbf P$$ and $$\mathbf M$$ both vanish, and also happens to contain the expression for the electromagnetic energy density $$u = \frac 1 2 \mathbf E \cdot \mathbf D + \frac 1 2 \mathbf H \cdot \mathbf B.$$ Maybe I'm missing some relevant vector calculus identities, or relations between the fields. Otherwise, I guess the only way out is to express $$\mathbf D$$ and $$\mathbf B$$ in terms of $$\mathbf E$$ and $$\mathbf H$$ respectively, through the usual relations (valid in all homogeneous isotropic linear media) $$\mathbf D = \varepsilon \mathbf E, \qquad \mathbf B = \mu \mathbf H,$$ with $$\varepsilon = \varepsilon_0 \varepsilon_r = \varepsilon_0 (1+\chi_e)$$ and $$\mu = \mu_0 \mu_r = \mu_0(1+\chi_m)$$, which would yield $$T_{ij} = \varepsilon E_i E_j + \frac 1 {\mu} B_i B_j - \frac 1 2 \delta_{ij} \left( \varepsilon E^2 + \frac 1 {\mu} B^2\right),$$ but this means we lose generality, and we won't be able to treat e.g. ferroelectric or ferromagnetic materials, or anisotropic linear media.

Any thoughts?