Given two Riemannian manifolds $M,N$, we say that $f:M \to N$ is harmonic if it is a critical point of the Dirichlet energy functional.
More precisely, this means that for every variation $f_t$ of $f$ with variation-field $v:=\frac{\partial f_t}{\partial t}|_{t=0}$ which is compactly supported in the interior of $M$, $\frac{d E(f_t)}{dt}|_{t=0}=0$.
Using the fact that this property of $f$ is equivalent to $f$ being a solution of a certain differential equation, one can deduce the following statement:
Claim: Suppose that for every point $p \in M$, there exist an open neighbourhood of $p$, $U_p\subseteq M$, such that $f|_{U_p}:U_p \to N$ is harmonic. Then $f$ is harmonic as a map $M \to N$.
Question: Is there a way to prove this claim without the passage $$\text{being critical} \to \text{satisfying E-L equation} \to \text{being critical}?$$
In other words, suppose you only know the "critical point definition" (and never heard of Euler-Lagrange equations). Is there a way to see directly that this property is local?
A naive idea is that given an arbitrary variation, we can somehow represent it as a finite number of compositions of "small variations" but I am not sure this makes any sense or really helpful.
For start, I am ready to assume $N=\mathbb{R}^n$ if it makes the problem easier.
Here is one way to solve this puzzle. Let ${\mathcal U}$ be the cover of $M$ by open subsets $U_p$ as in your question. Since $M$ is paracompact, without loss of generality, we may assume that this cover is locally finite. In particular, it admits a partition of unity $\{\eta_U\}_{U\in {\mathcal U}}$. The manifold $N$ admits an isometric embedding in some Euclidean space ${\mathbb R}^k$. Henceforth, I will consider $N$ as a submanifold of $E^k$ with the induced Riemannian metric. Then the energy of a map $$ f: M\to N $$ is the same as the energy of $f$ as a map $f: M\to {\mathbb R}^k$. One has to be careful here because $M$ is noncompact: We will have to take the energy of the restriction of $f$ to a relatively compact subset of $M$. To simplify the matters, let me simply assume that $M$ is compact, otherwise I have to repeatedly restrict to arbitrary relatively compact subsets.
For the purpose of variation of the energy functional we, of course, still have to restrict to maps $M\to N$. On the infinitesimal level, we will be computing $$ \delta E(f)(V) $$ where $V$ is a tangent vector field along $N$. It is convenient to regard $V$'s as defined on $M$ rather than on $N$, so technically speaking, they are sections of the pull-back of the tangent bundle: $f^*(TN)$. I will denote the space of such sections as ${\mathfrak X}_f(M,N)$.
Notice that I have no need to actually know the explicit form of the derivative of $E$. In particular, I do not even need to know what energy is, as long as $E$ is an integral $\int_M e(f)$, where $e(f)$ (which normally is the energy-density) is some (nonlinear) 1st order differential operator on maps $M\to {\mathbb R}^k$ (but I need to know that $E$ is differentiable in Gateaux's sense).
The condition that $f: M\to N$ is a critical point of $E$ is that for every tangent vector field $V$ along $N$ we have $$ \delta E(f)(V)=0. $$
Notice that $\delta E(f)(V)$ is linear in the variable $V$.
Now, suppose that I know only the vanishing property
$$ \delta E(f|_{U})(V)=0, \forall V\in {\mathfrak X}_f(U, N) $$
for each $U\in {\mathcal U}$. Given a vector field $V\in {\mathfrak X}_f(M, N)$, I define local vector fields $$ V_U:= \eta_U V $$ supported on $U$. Since $\{\eta_U\}$ is a partition of unity, we have $$ V= \sum_{U\in {\mathcal U}} V_U. $$ Then, by linearity, $$ \delta E(f)(V)= \sum_{U\in {\mathcal U}} \delta E(f|U)(V_U). $$ Each term in the right hand side of this equation vanishes since $f$ was a "local" critical point of $E$. Hence, $\delta E(f)(V)=0$ for every $V$ as above.