Implement pointwise observation operator

[finsberg]

Adding function to construct pointwise observation matrix $\mathbf{B}$ than maps a function $u$ to points in $\mathbb{R}^m$, i.e

$$ \mathbf{d} = \mathbf{B}\mathbf{u} $$

If $u$ is a function

$$ u(x) = \sum_{j=1}^N \mathbf{u}_j \phi_j(x) $$

then the entries are defined as:

$$ \mathbf{B}_{ij} = \phi_j(x_i) $$

This is useful if you e.g have a misfit functional $J$ given by

$$ J(\mathbf{u}) = \frac{1}{2} | \mathcal{O}(\mathbf{u}) - \mathbf{d}_{obs} |^2 $$

where $\mathcal{O}$ is an operator that maps $u$ to observations (for examples if we have observations of $u$ at specific points). If $\mathcal{O}(\mathbf{u}) = \mathbf{B}\mathbf{u}$, then the gradient is simply

$$ \nabla J(\mathbf{u}) = \mathbf{B}^T (\mathbf{B}\mathbf{u} - \mathbf{d}_{obs}) $$

[jorgensd]

@finsberg i think you could use: https://github.com/scientificcomputing/fenicsx_ii/blob/main/src/fenicsx_ii/interpolation.py where you make your points into a 0-D mesh with a DG-0 space.
Especially as the current code assumes an identity transformation, no piola mapped elements can be used.