In single dimension calculus, there is one and only one obvious way of integrating: over an interval. Of course, there are multiple ways of defining the integral (Riemann sum, Darboux integral, Lebesgue integral), but, where defined, these definitions agree.
Once we move from $\mathbb R$ to $\mathbb R^n$, there seems to be a plethora of types of integrals, with most sources making no attempt to unify or organize them:
In this section we will continue looking at line integrals and define the second kind of line integral we’ll be looking at... In this section we will define the third type of line integrals we’ll be looking at: line integrals of vector fields.
The Area under a Curve and its Many Generalizations... The method of doing this used is generalized to define a wide variety of integrals...
Is there a way to unify or at least organize all these types of integrals?
My attempt to do so is below. Is it correct? What revisions does it need? Update: My goal is to do this in a simple way, without invoking the complexity and abstraction of differential forms.
For $f: \mathbb R \to \mathbb R^m$, the one obvious way to integrate is over an interval of $\mathbb R$. When instead $f$ takes an input in $\mathbb R^n$, we gain one, and in some cases two, new ways to integrate:
Region ("Multiple") Integrals
For multiple dimensions, with $f: \mathbb R^n \to \mathbb R^m$, we can take the integral over any "nice" $k \leq n$ dimensional region of $\mathbb R^n$. We can do this for regions of any dimension $k \leq n$ by partitioning it into very small subregions. However, for $k > 1$, these subregions are not defined by two endpoints such that we can define $\Delta x_i$. Rather, we need to take the size or region measure (hypervolume, volume, area, etc.) appropriate for $k$ (irrespective of $n$) of each subregion. This partition is, unlike intervals, unordered and unoriented, and the size of each subregion is always positive: $$\sum_{r_{i, j,...} \in R}f(\xi_{i,j,...})\cdot \|r_{i,j,...}\|.$$
The integral is the limit of such as $\max_{i,j,...} \|r_{i,j,...}\| \to 0$ (provided it exists). (And the fact that the partition unordered and unoriented, with each subregion having positive size, is the reason why the Change of Variables Theorem uses the absolute value of the determinant and not the determinant itself.)
By Fubini's theorem, this limit is (under reasonable conditions) equal to an iterated integral, and so these are often called multiple integrals. However, it is important to emphasize that there is nothing "multiple" in the definition of multiple integrals: They are limits of single Riemann sums over partitions of regions with dimension $> 1$.
This definition therefore includes volume integrals, area integrals, line integrals over scalar fields; all are simply region integrals, defined the same way, with the region having different dimensions.
Vector Field Integrals: Line (dim = 1) and Flux (codim = 1)
Furthermore, if $n = m$, then $f$ defines a vector field, and we gain the ability to take the dot product $f(x) \cdot x$. This gives us two new ways to take integrals: line integrals and flux integrals.
If the region is a curve (i.e. $k = 1$), we can approximate each subregion $r_i$ by a vector $s_i \in \mathbb R^n$ and take the Riemann sum of $$\sum f(\xi_i) \cdot s_i.$$
And if the region is a surface (i.e. $k = n - 1$), we can construct a flux integral by approximating each $r_i$ by an $k$ dimensional hyperplane $p_i$ and letting $n_i$ be the unique vector that is normal to $p_i$, has magnitude equal to the "size" (volume, area, etc.) of $p_i$, and has a sign determined by convention, giving the Riemann sum of $$\sum f(\xi_i) \cdot n_i.$$
I quite like Terry Tao's comment/explanation that the distinct notions of 'integral'/'integration' are really and essentially about distinct things:
That they coalesce so nicely in the one-dimensional case may very well be a small miracle