Historically, it was Maxwell who introduced the objects and Faraday - I believe - who invoked the idea behind them of treating the electromagnetic field as a dynamic continuum pervading all of space. They are simply the objects that are central to converting the equations for power density $𝔓$ and force density $𝕱$ (which already pertain to continuua):
$$𝔓 = 𝐉·𝐄, \hspace 1em 𝕱 = ρ𝐄 + 𝐉×𝐁$$
fully to continuum form, by retro-applying the Maxwell equations
$$∇·𝐃 = ρ, \hspace 1em ∇×𝐇 - \frac{∂𝐃}{∂t} = 𝐉$$
for the charge density $ρ$ and current density $𝐉$ to replace them by the response fields $𝐃$ and $𝐇$.
This has the same utility as rewriting Newton's laws in continuum form as transport equations in fluid dynamics. It bears mentioning that historically, the continuum form came first even in Newton's laws. He expressed his laws, in The Principia introducing the concept of mass density first, not mass. So, this applies anywhere you need to do things in terms of continuua and fluid dynamics. That includes, as a special case, the setting of General Relativity, since Einstein's Field Equations are - essentially - fluid dynamics equations with law of gravity woven into them.
A relevant result, for instance, could be posed that if there is at least one frame of reference in which a constitutive law exists that relates the response fields to the force fields $𝐄$ and $𝐁$ through some generating function $𝔏(𝐁,𝐄)$ by:
$$𝐃 = \frac{∂𝔏}{∂𝐄}, \hspace 1em 𝐇 = -\frac{∂𝔏}{∂𝐁}$$
and if the Maxwell equations
$$∇·𝐁 = 0, \hspace 1em ∇×𝐄 + \frac{∂𝐁}{∂t} = 𝟬$$
hold for the force fields, then you can re-express the force and power densities (in that frame of reference) entirely in terms of various scalar, vector and dyad densities by:
$$\begin{align}
-𝔓 &= ∇·(𝐄×𝐇) + \frac{∂(𝐃·𝐄 - 𝔏)}{∂t}, \\
𝕱 &= ∇·(𝐃 𝐄 + 𝐁 𝐇) - ∇(𝐁·𝐇 + 𝔏) + \frac{∂(𝐁×𝐃)}{∂t}
\end{align}$$
using dyad notation to denote the terms $𝐃 𝐄 + 𝐁 𝐇$ in the stress tensor. So, now you have transport laws in the framework of fluid dynamics.
You can prove this by straightforward substitution, the use of the Maxwell equations for $(𝐁,𝐄)$, the use of the chain rule on $𝔏$ and some vector algebra.
This result holds, unaltered, in both the Relativistic and non-relativistic settings, because the four Maxwell equations, cited above, used in the derivation, reside at a deeper layer of the underlying space-time geometry where distinctions between space-like versus time-like, causal structure, relativitic versus non-relativistic, scale, orthogonality, etc. are not seen. It's "pre-metrical", to use a term Hehl has invoked in his writings.
Maxwell didn't actually use any derivation as simple as this, because he never got a proper handle on what the force and power laws ought to have been, though there clues were staring him in the face even at the time. Hint: consider how force and power density transform under Galilean boosts ... especially relevant since he made a whole point in treating and discussing the transform laws for his other equations. But this was the general result he was aiming for.
His analysis is recounted in chapter XI of the second book in his treatise, as seen over here
ON ENERGY AND STRESS IN THE ELECTROMAGNETIC FIELD.
https://en.wikisource.org/wiki/A_Treatise_on_Electricity_and_Magnetism/Part_IV/Chapter_XI
(Note: his $μ$ and $𝐇$ differ from our $μ$ and $𝐇$ by factors of $4π$ and his $𝐄$ differs from our $𝐄$ by the inclusion of a velocity-dependent term ... which has to do with the "stationary" versus "moving" forms of the constitutive laws, as alluded to below.)
The special case, that you're already more familiar with, and consistent with his analysis, is where the generating function is of the form:
$$𝔏(𝐁,𝐄) = \frac{ε E^2} 2 - \frac{B^2}{2μ},$$
which gives rise to the constitutive relations $𝐃 = ε 𝐄$ and $𝐁 = μ 𝐇$. In the non-relativistic setting, and in the general setting of media other than a relativistic vacuum, these equations generally hold in only one frame of reference, and so - historically - they (when combined with the rest of the Maxwell equations) were referred to as the "stationary form" of the equations; the "moving form" would require the inclusion of a velocity $𝐆$. Maxwell's $𝐄$, in this setting, would be equivalent to our $𝐄 + 𝐆×𝐁$. The "moving" form of the constitutive laws are
$$𝐃 + α 𝐆×𝐇 = ε(𝐄 + 𝐆×𝐁), \hspace 1em 𝐁 - α 𝐆×𝐄 = μ(𝐇 - 𝐆×𝐃)$$
and are the (non-relativistic) relations of Maxwell and Thomson if $α = 0$, and are equivalent to the equations Lorentz posed for moving media in the early 1900's. The $- 𝐆×𝐃$ was Thomson's correction and was not present in Maxwell's analyses.
They are the (Relativistic) Maxwell-Minkowski relations, derived by Minkowski in his seminal 1908 "Minkowski geometry" paper and independently by Einstein and Laub in 1908, if $α = (1/c)^2$.
If $εμ = α$, the "moving" forms are (almost) equivalent to the "stationary" forms. In particular, for Relativity, in a vacuum, where $ε = ε_0$ and $μ = μ_0$, with the wave speed $1/\sqrt{εμ} = 1/\sqrt{ε_0μ_0} = c$, the constitutive laws can be asserted in all inertial frames - both "stationary" and "moving". (That was Einstein's key point in his 1905 paper on the electrodynamics of moving bodies).
It's actually a bit easier to use the Routhian
$$ℜ(𝐄,𝐇) = 𝐁·𝐇 + 𝔏(𝐁,𝐄)$$
which assumes that the Hessian dyad $∂^2𝔏/∂𝐁∂𝐁$ is non-singular, since it's more consistent with Maxwell's terming of $𝐁$ and $𝐃$ as somehow "induced" from the "force" fields $𝐄$ and $𝐇$. In that case, the constitutive law would become
$$𝐃 = \frac{∂ℜ}{∂𝐄}, \hspace 1em 𝐁 = \frac{∂ℜ}{∂𝐇},$$
and the relevant result would be expressed as
$$-𝔓 = ∇·(𝐄×𝐇) + \frac{∂(𝐃·𝐄 + 𝐁·𝐇 - ℜ)}{∂t}, \hspace 1em 𝕱 = ∇·(𝐃 𝐄 + 𝐁 𝐇) - ∇ℜ + \frac{∂(𝐁×𝐃)}{∂t},$$
and for the previously-mentioned special case, we'd have:
$$ℜ(𝐄,𝐇) = \frac{ε E^2} 2 + \frac{μ H^2} 2.$$
This is closer to the expressions Maxwell used in his analyses.
If you want to do the analysis correctly and more comprehensively, you should also be accounting for the source terms. J. D. Jackson did that in his classical book on Electromagnetism. Here, the generating function
$𝔏(φ,𝐀,𝐁,𝐄)$ is now assumed to be not just a function of the force fields $(𝐁,𝐄)$, but also of the potentials $(φ,𝐀)$ which they are derived from by the Maxwell equations:
$$𝐁 = ∇×𝐀, \hspace 1em 𝐄 = -∇φ - \frac{∂𝐀}{∂t}.$$
The constitutive laws are then expanded to also account for the charge and current densities by:
$$ρ = -\frac{∂𝔏}{∂φ}, \hspace 1em 𝐉 = \frac{∂𝔏}{∂𝐀}.$$
Then, with the aid of the equation of continuity,
$$∇·𝐉 + \frac{∂ρ}{∂t} = 0$$
which we have to assume hold for these sources, for consistency, both the left and right-hand sides of the force and power laws can be pulled into continuum form and expressed as equations of continuity in their own right:
$$\begin{align}
0 &= ∇·(𝐄×𝐇 - 𝐉φ) + \frac{∂(𝐃·𝐄 - 𝔏 - ρφ)}{∂t},\\
𝟬 &= ∇·(𝐃𝐄 + 𝐁𝐇 + 𝐉𝐀) - ∇(𝐁·𝐇 + 𝔏) + \frac{∂(𝐁×𝐃 + ρ𝐀)}{∂t}.
\end{align}$$
