Are all symmetries of the Dirichlet functional isometries? This is a cross-post from MSE (no answer there).
Let $M,N$ be oriented $d$-dimensional Riemannian manifolds, $M$ compact*, and let $f:M \to N$ be smooth.
Consider the Dirichlet  energy functional: $E_{M,N}(f)=\int_M \|df\|^2 \operatorname{Vol}_g$.
($\operatorname{Vol}_g$ is the Riemannian volume form of $M$).
Fix another $d$-dimensional Riemannian manifold $W$.
For every isometric immersion $\phi:N \to W$,
$$ E_{M,N}(f)=E_{M,W}(\phi \circ f) \, \text{ for any } \, f:M \to N. \tag{1}$$
Let us call a smooth map $\phi$ which sarisfies $(1)$ a  symmetry of the Dirichlet's integral.

Question: Is every symmetry an isometric immersion? 
Does anything change if we restrict the symmetries to be invertible maps? (or diffeomorphisms)? 

Has this notion of symmetry been studied in the context of mappings between manifolds? (In more "classical Euclidean" settings I found similar notions by the name of "variational symmetry groups").
*We can relax the compactness assumption, but then we need to restrict the discussion to maps which are constant outside a compact domain.

Comments:
$(1)$ The "pointwise analog" is trivial:
If $B \in M_d$, and for every $A \in M_d$, $ \|BA\|=\|A\|$ (where $\| \|$ is the Euclidean norm), then $B$ is an orthogonal matrix.
The challenge is that the notion of a Dirichlet's-symmetry is a global one, while the concept of isometric immersion is local. 
I think the rough idea should be to choose "test maps" $f$ which are very "localized" (are constant outside small balls). However, this does not seem trivial, since the differential would have to pass from a given linear map to zero, so its norm would vary. (Think of a bump function which goes quickly from one to zero, you can't make the integral of the derivative small).
$(2)$ The answer could a-priori depend on the manifolds. Even the case where $M$ is an Euclidean ball in $\mathbb{R}^d$, $N=W=\mathbb{R}^d$, the answer does not seem to be trivial (see the previous comment).
 A: The answer to your question is "yes, every symmetry (in the sense you have specified) is an isometric immersion".
To see why, first note that, if $(M,g)$ is a compact Riemannian $m$-manifold and $h$ is any smooth quadratic form on an $n$-manifold $N$, one can always define the Dirichlet energy $E_{g,h}(f)$ of a smooth map $f:M\to N$.  You don't need that $h$ be positive definite (or even nondegenerate) to do this.  (You don't need $M$ to be oriented either; instead of using the $g$-volume, you can use the $g$-density instead.)  To see this, just look at the formula for the Energy density in local coordinates $x=(x^i)$ on $M$ and $y = (y^a)$ on $N$.  In those coordinates, you'll have $g = g_{ij}(x)\,\mathrm{d}x^i\mathrm{d}x^j$ and $h = h_{ab}(y)\,\mathrm{d}y^a\mathrm{d}y^b$ and, when the map is expressed as $y = f(x)$, the Dirichlet energy density becomes
$$
\Phi_{g,h}(f) = h_{ab}\bigl(f(x)\bigr)\frac{\partial^af}{\partial x^i}\frac{\partial^bf}{\partial x^j} g^{ij}(x)\,\,G(x)\,\mathrm{d}x\tag 1
$$
where $G^2 = \det(g_{ij})$.  Note that, while you do need $g$ to be nondegenerate, you don't need $h$ to be nondegenerate.  (I'm omitting the traditional coefficient of $\tfrac12$ since you did.)
As a result, your question reduces to the following one:  Suppose that $E_{g,h}(f) = 0$ for all $f:M\to N$.  Does this imply that $h$ vanishes identically?  
Now, this is a question of the form of so-called 'null Lagrangians', i.e., those first-order Lagrangians for which the Euler-Lagrange equations are trivial.  From the above hypothesis, $\Phi_{g,h}$ would be a null Lagrangian.  
It is a classical fact (see any decent book on the calculus of variations) that a null Lagrangian of first order from an $m$-dimensional manifold to an $n$-dimensional manifold must be of the local form
$$
\Phi(f) = \sum_{I,A} h^I_A\bigl(x,f(x)\bigr)D^A_I\left(\frac{\partial f}{\partial x}\right)\,\,\mathrm{d}x,\tag 2
$$ 
where the sum is over all multi-indices $I = (i_1,\ldots,i_k)$, $A=(a_1,\ldots,a_k)$
and $0\le k\le \min(m,n)$, where $1\le i_1 < i_2< \cdots < i_k\le m$ and $1\le a_1 < a_2<\cdots<a_k\le n$ and $D^A_I(L)$ is the corresponding $k$-by-$k$ minor of the $n$-by-$m$ matrix $L$.
However, one immediately sees that the above Dirichlet energy (1) is of the form (2) if and only if $h_{ab}(y)\equiv0$, which implies the desired result. 
A: This is just an elaboration on Robert's great answer:
The key idea is to use the fact that the induced metric on the Hom -space of two inner product spaces is "linear" in the "metric" on the target.
(The "metric" here means to any quadratic form, it does not need to be positive or even non-degenerate).
Denote our manifolds by
$(M,g),(N,h),W,\eta)$. Let $\phi$ be a (smooth) symmetry. Then $\phi:(N,\phi^*\eta) \to (W,\eta)$ is "isometric" (preserves the quadratic form), so $$E_{(M,g),(N,\phi^*\eta)}(f)=E_{(M,g),(W,\eta)}(\phi \circ f) \,  \, \text{ for any } \, f:M \to N.  \tag{1}$$
Combining this with
$$E_{(M,g),(N,h)}(f)=E_{(M,g),(W,\eta)}(\phi \circ f) \, \text{ for any } \, f:M \to N, \tag{2}$$
we obtain
$$ E_{(M,g),(N,h)}(f)=E_{(M,g),(N,\phi^*\eta)}(f),$$
or (using the said "linearity" ) $$ E_{(M,g),(N,h-\phi^*\eta)}(f)=0 \, \text{ for any } \, f:M \to N. \tag{3}$$
Hence, by Robert's argument $h-\phi^*\eta=0$, so $\phi$ is an isometric immersion.
