An action integral
$$
S = \int_{\Omega}L(x_1,...,x_n,f,D_{x_1}f, ...,D_{x_n}f) dx_1\cdots dx_n,
$$
can be minimized over a compact set $\Omega$ via the multi-dimensional generalization of Euler-Lagrange equation in the form
$$
\dfrac{\partial L}{\partial f} - \sum_{i=1}^{n} \dfrac{\partial }{\partial x_i} \dfrac{\partial L}{\partial D_{x_i}f} = 0
$$
with $D_{x_i}f := \partial f/\partial x_i$, $f : \mathbb{R}^n\to \mathbb{R}$, and Lagrangian $L(x_1,...,x_n,f,D_{x_1}f, ...,D_{x_n}f)$.
I still couldn't figure out the derivation/proof of this multi-d extension. First of all, could someone point me to a nice derivation/proof of this extended case?
Having said that, I am particularly wondering how we could adopt this equation for arbitrary domains beyond product topology of intervals (i.e., $\Omega = I_1 \times I_2 \times \cdots \times I_n$ for closed intervals $I_i \in \mathbb{R}$). For example, when the domain is a closed n-ball $\Omega = B^n$. Do the same equations hold for any sufficient set of boundary conditions defined on the boundary of the ball $\partial \Omega = S^n$?
The equation you wrote is the generalization of the usual Euler-Lagrange equation from classical mechanics to classical field theory. You can find the derivation of this in a lot of places, just try searching for the "classical field theory version" of the Euler-Lagrange equations. For instance, you can look at Section 5.3 of Nastase's Classical field theory.
Also, there's no reason that these equations can only apply to products of intervals; the derivation doesn't assume anything like $\Omega = I_1 \times \cdots \times I_n$. Indeed, these equations hold for any compact manifold $\Omega$ with boundary (for instance, $\Omega = \mathsf{B}^n \subset \mathbb{R}^n$), as long as you place appropriate boundary conditions on the fields (i.e., the functions $f$) that you're considering. This is because the derivation involves some kind of "integration by parts" that introduces a boundary term, which you'd like to make zero in order to obtain the equations above.
Edit. To be clear, maybe you're under the impression that we can only work with products $\Omega = I_1 \times \cdots \times I_n$ since for the $1$-dimensional Euler-Lagrange equations (i.e. where we replace $x_1,\dots,x_n$ by $t$), we work with a single interval $\Omega = I$. But this is just because any compact $1$-dimensional manifold (connected, with nonempty boundary) has to be (diffeomorphic to) a closed interval! However, in dimensions $n > 1$, there are way more compact manifolds (with boundary) than just products of intervals.