Approximate Fréchet derivative of PDE solution operator Let's assume that $F$ is the solution operator of some PDE problem, i.e. $F$ maps some parameter or boundary condition function to the solution of a given PDE. Let also be $\Pi$ be the evaluation operator in a given point. How do I obtain an expression (or a numerical scheme for computation) of the Fréchet derivative of $\Pi\circ F$ in some point?
Or very concretely, consider the linear elliptic problem
$$-\nabla \cdot (k(x)\cdot \nabla p(x)) = 0 $$
with boundary conditions $p  = h$ on the boundary.
$F$ maps $k\mapsto p$ and $\Pi$ evaluates at $x_0$, i.e. $\Pi u = u(x_0)$. Then $(\Pi\circ F)(k) = p(x_0)$ with $p$ solving the PDE with parameter $k$.
What is the Fréchet derivative $D(\Pi\circ F)(k)$?
 A: Since point evaluation and the trace operator are linear, you only need to worry about the solution mapping $F:k\mapsto p$; the rest follows from the chain rule for Fréchet derivatives. In general, and for this PDE in particular, Fréchet differentiability of the solution mapping is a quite subtle issue, and depends on the exact function spaces your mapping is defined on.
But you are only asking about the representation of the derivative, not its existence, so we can just proceed formally. Your PDE defines an implicit function $p$ of $k$ via $e(p,k)=0$, so the derivative can be computed using the implicit function theorem: If $e(p,k)=0$, then so is the total derivative
$$ e_p(p,k)p' + e_k(p,k)k'=0,$$
where $e_p(p,k)$ denotes the partial derivative of $e$ with respect to $p$ etc.
In your concrete setting, this translates to
$$ -\nabla\cdot(k \cdot\nabla p') -\nabla\cdot(k' \cdot\nabla p) = 0.$$
Therefore, the derivative $p':=F'(k)k'$ acting on a direction $k'$ is given by the solution of
$$\left\{\begin{aligned} -\nabla\cdot(k \cdot\nabla p') &= \nabla\cdot(k' \cdot\nabla p) &&\text{on }\Omega,\\ p' &= 0 &&\text{on }\partial\Omega,
\end{aligned}\right.$$
where $p$ is the solution of the original PDE for the parameter $k$. (The boundary conditions have to be homogeneous since the derivative is by definition a linear operator, and not affine.) This is a linear PDE that can be solved by standard methods, e.g., finite elements. The derivative of the composite mapping $\Pi\cdot F$ is then just the corresponding point evaluation of the trace of $p'$.
The next step is to insert this representation in the definition of the Fréchet-derivative 
$$ \lim_{\|k'\|\to 0} \frac{\|F(k+k')-F(k)-F'(k)k'\|}{\|k'\|}$$
and use suitable a priori estimates for the original and the linearized PDE to show convergence to zero, which of course only works for appropriate choices of the norms in the numerator and denominator. In this case, you can show differentiability from $L^\infty(\Omega)$ to $H^1(\Omega)$; see, e.g., 
Lemma 3.1 in
Ute Aßmann and Arnd Rösch, MR 3042664 Identification of an unknown parameter function in the main part of an elliptic partial differential equation, Z. Anal. Anwend. 32 (2013), no. 2, 163--178.
(An alternative is to verify that the conditions of a "proper" implicit function theorem in function spaces are satisfied.)
The general approach using the implicit function theorem to obtain a representation of the derivative is discussed in Section 1.6 of
M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich, MR 2516528 Optimization with PDE constraints, ISBN: 978-1-4020-8838-4 
or Section 3.2 of 
Juan Carlos De los Reyes, MR 3308473 Numerical PDE-constrained optimization, Springer Briefs in Optimization ISBN: 978-3-319-13394-2; 978-3-319-13395-9.
