Inverting differential equation using finite element methods

Question

tl;dr; How to use FEM tools to extract models needed to invert PDEs.

Context

In astrophysics, one is interested in so-called 'cosmic archeology' which involves recovering the origin of a given observation, while modelling its evolution. The idea is to be able to understand what may have caused in the past an given signature found in the data. For instance, can we explain the observed vertical velocity distribution of stars above and below the galactic disc seen by the Gaia spacecraft in terms of past satellites having hit our Milky way.

Example

As a test example let us consider a 1D diffusion equation sourced by a finite set of heat sources

 source[x_, t_] = 
 3 BSplineBasis[3, t 4] BSplineBasis[3, (x - 1/8) 4] +
   2 BSplineBasis[3, -2 + t 4] BSplineBasis[3, (x - 5/8) 4] +
   BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 2/8) 4] +
   BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 1/8) 4] +
   BSplineBasis[3, -1/2 + t 4] BSplineBasis[3, (x - 4/8) 4] +
   3/2 BSplineBasis[3, -3 + t 4] BSplineBasis[3, (x - 1/8) 4];

ContourPlot[source[x, t], {x, 0, 1}, {t, 0, 1}, PlotPoints -> 36, 
 Exclusions -> None, PlotRange -> All, 
 AspectRatio -> 1, Contours -> 10]

The diffusion diagram will obey

    sol0 = NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] == 
        source[x, t],
       f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0}, f, {x, 0, 1}, {t, 0, 2}];

     ContourPlot[sol0[x, t], {x, 0, 1}, {t, 0, 1}, FrameLabel -> {x, t}, 
 AspectRatio -> 1, PlotRange -> All, Contours -> 30, PlotPoints -> 50]

Here I have assumed arbitrarily that the edges of [0,1] did not let heat diffuse. I also assumed that initially there was no heat.

Let me first generate the corresponding data set of positing and time $(x,t)$ for later use

data = N[Flatten[
Table[{x, t, sol0[x, t]}, {x, 0, 1, 1/32}, {t, 0, 1, 1/32}], 1]];

My purpose is to invert this data set to recover the source of heating.

I other words, can I recover the first plot from the second one, if I assume I know how the heat source diffuse?

Attempt

I can define a set of spline functions which cover the $(x,t)$ space as follow:

nn = 16;
knots = Flatten[{{0, 0}, (Range[0, nn]/nn), {1, 1}}];

basis0 = Flatten@
   Table[BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, 
      t], {i, 0, nn}, {j, 0, nn}];

For instance, the 60th function obeys:

Plot3D[basis0[[60]], {x, 0, 1}, {t, 0, 1}, PlotRange -> All]

The image of this basis satisfies

basis = Flatten@
   Table[NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] == 
       BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, t],
      f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0}, 
     f[x, t], {x, 0, 1}, {t, 0, 1}], {i, 0, nn}, {j, 0, nn}];

Plot3D[basis[[60]], {x, 0, 1}, {t, 0, 1}, PlotRange -> All]

From this basis and the data I can generate the linear model a0 and a relating resp. the spline coefficients to the source map and the spline coefficients heat map:

ff = Function[{x, t}, basis0 // Evaluate];
a0 = ff @@ # & /@ (Most /@ data0);

and

ff = Function[{x, t}, basis // Evaluate];
a = ff @@ # & /@ (Most /@ data);
a // Image // ImageAdjust

Let me first check that I can fit the source map with the splines:

fit0[x_, t_] = 
basis0.LinearSolve[Transpose[a0].a0, Transpose[a0].(Last /@ data0)];
ContourPlot[fit0[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20, 
   PlotRange -> All]

Similarly, I can define an (isotropic) penalty corresponding to $\int |\Delta T|^2 dx dt$ as

ff = Function[{x, t}, D[basis0, x, x] + D[basis0, t, t] // Evaluate];
s0 = ff @@ # & /@ (Most /@ data0);
pen = Transpose[s0].s0; pen /= Max[Flatten[Abs[pen]]];
pen // Image // ImageAdjust

A solution to the inverse problem then follows simply from inverting a with a small roughness penalty as

sol[x_, t_] = 
  basis0.LinearSolve[Transpose[a].a + 10^-7 pen, 
    Transpose[a].(Last /@ data)];

ContourPlot[sol[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20, 
 PlotRange -> All]

Question

I am fairly certain my present implementation is effectively redundant with the way NDSolve can actually solves the differential equation using Finite Element methods. Hence my question:

Could one make use of the actual solver in NDSolve to formulate the inverse problem?

In other words, can we extract from the FEM toolkit FEM_a,FEM_source and FEM_solution and FEM_grid so that

   FEM_solution =  FEM_a  FEM_source

where 'FEM_' stands for as sampled by the underlying mesh,FEM_grid of the FEM toolkit.

This would be of interest in terms of efficiency, but also in order to address more complex and realistic inverse problems?

For Gaia data, the diffusion is in fact occurring in 3D and is anisotropic, so a robust and efficient formulation would help!

Technically I believe FEM have access to both a and pen so it would be great to be able to access them for the sake of solving the inverse problem.

I am guessing that this link would be a good starting point?

Comment

Note that the above implementation is partially incorrect at the top edge, because most spline basis elements are required to be zero on the boundary, whereas the correct solution should have an outgoing flux condition. This is something the FEM would handle naturally when the boundary conditions are taken care of.

  Plot[{sol[x, 1], sol0[x, 1]}, {x, 0, 1}]

Note that in astronomy we unfortunately don't have access to the full diffusion diagram but only typically to a given snapshot (i.e. data on a line at fixed time, and/or possibly the time derivative on that line). So we can only extrapolate in the past up to some quite limited time horizon.

Complement 1: 1+1D code

source[x_, t_] = 
  3 BSplineBasis[3, t 4] BSplineBasis[3, (x - 1/8) 4] +
   2 BSplineBasis[3, -2 + t 4] BSplineBasis[3, (x - 5/8) 4] +
   BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 2/8) 4] +
   BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 1/8) 4] +
   BSplineBasis[3, -1/2 + t 4] BSplineBasis[3, (x - 4/8) 4] +
   3/2 BSplineBasis[3, -3 + t 4] BSplineBasis[3, (x - 1/8) 4];
sol0 = NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] == 
     source[x, t],
    f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0}, 
   f, {x, 0, 1}, {t, 0, 2}];
nn = 16; knots = Flatten[{{0, 0}, (Range[0, nn]/nn), {1, 1}}];
basis0 = Flatten@
   Table[BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, 
      t], {i, 0, nn}, {j, 0, nn}];
basis = Flatten@
   Table[NDSolveValue[{D[f[x, t], t] - 1/4 D[f[x, t], x, x] == 
       BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, t],
      f[x, 0] == 0, f[0, t] == 0, f[1, t] == 0}, 
     f[x, t], {x, 0, 1}, {t, 0, 1}], {i, 0, nn}, {j, 0, nn}];
data = N[Flatten[
    Table[{x, t, sol0[x, t]}, {x, 0, 1, 1/32}, {t, 0, 1, 1/32}], 1]];
data0 = N[
   Flatten[Table[{x, t, source[x, t]}, {x, 0, 1, 1/32}, {t, 0, 1, 
      1/32}], 1]];
ff = Function[{x, t}, basis0 // Evaluate];
a0 = ff @@ # & /@ (Most /@ data0);
ff = Function[{x, t}, basis // Evaluate];
a = ff @@ # & /@ (Most /@ data);
fit0[x_, t_] = 
  basis0.LinearSolve[Transpose[a0].a0, 
    Transpose[a0].(Last /@ data0)];
fit[x_, t_] = 
  basis.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)];
ff = Function[{x, t}, D[basis0, x, x] + D[basis0, t, t] // Evaluate];
s0 = ff @@ # & /@ (Most /@ data0);
pen = Transpose[s0].s0; pen /= Max[Flatten[Abs[pen]]];
sol[x_, t_] = 
  basis0.LinearSolve[Transpose[a].a + 10^-7 pen, 
    Transpose[a].(Last /@ data)];
ContourPlot[source[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20, 
 PlotRange -> All,Exclusions -> None]
ContourPlot[sol[x, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20, 
 PlotRange -> All]

Complement 2: 2+1D codes

For the sake of completeness and to demonstrate why a more efficient implementation is needed here is the code for 2D diffusion without FEM (which for n=16 would take a white to run!).

source[x_, y_, t_] = BSplineBasis[3, t ] BSplineBasis[3, x]*
  BSplineBasis[3, y]  
sol0 = NDSolveValue[{D[f[x, y, t], t] - 1/4 D[f[x, y, t], x, x] - 
     1/4 D[f[x, y, t], y, y] == source[x, y, t], f[x, y, 0] == 0,
   DirichletCondition[f[x, y, t] == 0, True]}, f, {x, 0, 1}, {y, 0, 1}, {t, 0, 1}]
nn = 2;knots = Flatten[{{0, 0}, (Range[0, nn]/nn), {1, 1}}];
basis0 = Flatten@
   Table[BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, y]
     BSplineBasis[{3, knots}, k, t], {i, 0, nn}, {j, 0, nn}, {k, 0, nn}];
basis = Flatten@(Table[
      ParallelTable[
       NDSolveValue[{D[f[x, y, t], t] - 1/4 D[f[x, y, t], x, x] - 
           1/4 D[f[x, y, t], y, y] == 
          BSplineBasis[{3, knots}, i, x] BSplineBasis[{3, knots}, j, 
            y] BSplineBasis[{3, knots}, k, t], f[x, y, 0] == 0,
         DirichletCondition[f[x, y, t] == 0, True]}, 
        f[x, y, t], {x, 0, 1}, {y, 0, 1}, {t, 0, 1}], {j, 0, nn}, {k, 
        0, nn}], {i, 0, nn}]);
  data0 = N[Flatten[Table[{x, y, t, source[x, y, t]}, {x, 0, 1, 1/nn/2}, 
           {y, 0,1, 1/nn/2}, {t, 0, 1, 1/nn/2}], 2]];
data = N[Flatten[
    Table[{x, y, t, sol0[x, y, t]}, {x, 0, 1, 1/nn/2}, {y, 0, 1, 
      1/nn/2}, {t, 0, 1, 1/nn/2}], 2]];
ff = Function[{x, y, t}, basis // Evaluate];
a = ParallelMap[ff @@ # &, Most /@ data];
ff = Function[{x, y, t}, D[basis0, x, x] +
 D[basis0, y, y] + D[basis0, t, t] // Evaluate];
s0 = ff @@ # & /@ (Most /@ data0);
pen = Transpose[s0].s0; pen /= Max[Flatten[Abs[pen]]];
sol[x_, y_, t_] = 
  basis0.LinearSolve[Transpose[a].a + 10^-9 pen, 
    Transpose[a].(Last /@ data)];
ContourPlot[sol[x, 1/2, t], {x, 0, 1}, {t, 0, 1}, Contours -> 20, 
 PlotRange -> All]

Complement 3: Background

Let

$$\mathcal{L}\psi = \rho $$

represent a (linear) partial differential equation (possibly time dependant). I will assume that there exist a basis function over which I can project $\psi$, so that $$ \psi(x)=\sum_n a_n \phi_n(x)\,,$$ where I also request that all $\phi_n(x)$ satisfy the boundary conditions of the partial differential equation ( $x$ can represent say $\mathbf{r}$ or $(\mathbf{r},t)$), i.e. the analysis is not necessary limited to stationary PDE). If I put this expansion into the PDE, multiply by $\phi_p(x)$ (or a Dirac function as a variant, see below) and integrate over $x$, I get formally $$ \mathbf{L}\cdot \mathbf{\Phi} = \mathbf{P}\,, $$ where $L_{ij}= \int d x \phi_i \mathcal{L} \phi_j $, $P_{i}= \int d x \phi_i \rho $ and ${\Phi}_i= a_i$.

I can then invert for $ \mathbf{\Phi}$ as $$ \mathbf{\Phi} =\mathbf{L}^{(-1)} \cdot\mathbf{P}\,, $$ where $\mathbf{L}^{(-1)}$ is the (possibly regularised) pseudo inverse of $\mathbf L$ (e.g. through least square). This is a possible method for solving PDEs. I am assuming (wrongly?) that linear FEM methods are a variant of this technique.

Conversely, If I start with the solved equation

$$\psi = \mathcal{L}^{-1}\rho \,. $$ I can expand $\rho$ over a basis function,$\rho=\sum_n a_n \rho_n$ , project as previously and write eventually

$$ \mathbf{P} =\mathbf{R}^{(-1)}\cdot \mathbf{\Phi}\,, $$ where $\mathbf{R}^{(-1)}$ is the (possibly regularised) pseudo inverse of $\mathbf R$, whose components are $R_{ij}= \int d x \rho_i \mathcal{L}^{-1} \rho_j $. In my code above I have implemented something closely related to the second method, relying on NDSolve (I use a Dirac function instead of $ \rho_i$ to simply sample the measured $\phi(x)$). The reason is I am after recovering the source $\mathbf P$ given some knowledge of the response $\mathbf \Phi$.

My hope is that since FEM method solve the first problem they should have in store the tools to solve the second problem more efficiently?

Answer 1

I am going to show how to write this part of your post

source[t_, x_] = 
  3 BSplineBasis[3, t 4] BSplineBasis[3, (x - 1/8) 4] + 
   2 BSplineBasis[3, -2 + t 4] BSplineBasis[3, (x - 5/8) 4] + 
   BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 2/8) 4] + 
   BSplineBasis[3, -1 + t 4] BSplineBasis[3, (x - 1/8) 4] + 
   BSplineBasis[3, -1/2 + t 4] BSplineBasis[3, (x - 4/8) 4] + 
   3/2 BSplineBasis[3, -3 + t 4] BSplineBasis[3, (x - 1/8) 4];

tEnd = 2;
AbsoluteTiming[
 sol0 = NDSolveValue[{D[f[t, x], t] - 1/4 D[f[t, x], x, x] == 
      source[t, x], f[0, x] == 0, f[t, 0] == 0, f[t, 1] == 0}, 
    f, {x, 0, 1}, {t, 0, tEnd}
    , Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"FiniteElement"}}
    ];]
(* {0.337181, Null} *)

with the low level FEM functions. It's still not quite clear to me how you want to make use of this. More on this later. Note that I added a method option to tell NDSolve to actually make use of the FEM. Not all of the NDSolve calls you show actually make use of the FEM. But I think the method used is also not relevant.

To understand the code that follows I'd advise to read the FEMProgramming tutorial.

Set up the region, bcs, ics:

region = Line[{{0}, {1}}];
bcs = {DirichletCondition[u[t, x] == 0, True]};
initialConditionValue = 0.;
vd = NDSolve`VariableData[{"DependentVariables" -> {u}, 
    "Space" -> {x}, "Time" -> t}];

Needs["NDSolve`FEM`"]
nr = ToNumericalRegion[region];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];

Set up the PDE coefficients without the load term:

dim = RegionDimension[region];
initCoeffs = 
  InitializePDECoefficients[vd, 
   sd, {"DampingCoefficients" -> {{1}}, 
    "DiffusionCoefficients" -> {{-1/4 IdentityMatrix[dim]}}}];

We omit the load term for now, as that is the term that is variable in your examples and we will take care of that later.

Initialize the BCs, method data and compute the stationary (time independent) discretization and boundary conditions of the PDE (without the load). These coefficients and discretizations are the same for all the PDEs you solve so we compute them only once.

initBCs = InitializeBoundaryConditions[vd, sd, {bcs}];
methodData = InitializePDEMethodData[vd, sd];

sdpde = DiscretizePDE[initCoeffs, methodData, sd, "Stationary"];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

Now, we want to write a residual function for NDSolve to time integrate. At the same time we want the source to be variable.

makeResidualFunction[load_] := With[
  {loadCoeffs = 
    InitializePDECoefficients[vd, 
     sd, {"LoadCoefficients" -> {{load}}}]},
  With[
   {sloaddpde = 
     DiscretizePDE[loadCoeffs, methodData, sd, "Stationary"]},
   discretizePDEResidual[t_?NumericQ, u_?VectorQ, dudt_?VectorQ] := 
    Module[{l, s, d, m, tloaddpde},

     NDSolve`SetSolutionDataComponent[sd, "Time", t];
     NDSolve`SetSolutionDataComponent[sd, "DependentVariables", u];

     {l, s, d, m} = sdpde["SystemMatrices"];

     (* discretize and add the laod *)
     (*l+=sloaddpde["LoadVector"];*)
     tloaddpde = 
      DiscretizePDE[loadCoeffs, methodData, sd, "Transient"];
     l += tloaddpde["LoadVector"];

     DeployBoundaryConditions[{l, s, d}, sbcs];

     d.dudt + s.u - l
     ]
   ]
  ]

This functions get the 'source' function and defines a residual function.

Generate the initial conditions with boundary conditions deployed.

init = Table[
   initialConditionValue, {methodData["DegreesOfFreedom"]}];
init[[sbcs["DirichletRows"]]] = Flatten[sbcs["DirichletValues"]];

Get the sparsity pattern for the damping matrix for the time integration.

sparsity = sdpde["DampingMatrix"]["PatternArray"];

Set up the residual function.

makeResidualFunction[source[t, x]]

Do the time integration

AbsoluteTiming[
 ufun = NDSolveValue[{
    discretizePDEResidual[t, u[t], u'[ t]] == 0
    , u[0] == init}, u, {t, 0, tEnd}
   , Method -> {"EquationSimplification" -> "Residual"}
   , Jacobian -> {Automatic, Sparse -> sparsity}
   (*,EvaluationMonitor\[RuleDelayed](monitor=Row[{"t = ",CForm[t]}])*)
   , AccuracyGoal -> $MachinePrecision/4, 
   PrecisionGoal -> $MachinePrecision/4
   ]
 ]

(* {0.429631,.... *)

As you see the time integration is somewhat slower from top level code.

Convert the result to an interpolating function:

ffun = ElementMeshInterpolation[{ufun["Coordinates"][[1]], 
   methodData["ElementMesh"]}, Partition[ufun["ValuesOnGrid"], 1]]

Check that this is on the same order as the NDSolve result.

Plot3D[sol0[t, x] - ffun[t, x], {t, 0, tEnd}, {x, 0, 1}, 
 PlotRange -> All]

Discussion:

I think you make an implicit assumption that is not correct. You assume that the matrix assembly process is the expensive part. But it's not. It's the actual time integration that you need to do many many times that is expensive. Precomputing the system matrices can probably save a little when the parallel computation is run but you can not make the time integration go away altogether.

Answer 2

Let me try and begin to answer my own question, as it might inspire better answers. Here I will solve the Poisson equation as a test case using 0-splines.

Needs["NDSolve`FEM`"];
reg0 = Rectangle[{0, 0}, {1, 1}];
mesh0 = ToElementMesh[reg0, MaxCellMeasure -> 0.025, AccuracyGoal -> 1]

I can then extract the mesh elements

idx = mesh0["MeshElements"][[1, 1]];mesh0["Wireframe"]

In order to define the density on each cell I need to find the convex hull of each cell

pol = Table[mesh0["Coordinates"][[ idx[[i]]]] // ConvexHullMesh, {i,Length[idx]}];

I can then use the function RegionMember to define the Indicator of that cell (as shown in this answer)

basis = Table[f[x_, y_] = Boole[ RegionMember[pol[[i]], {x, y}]]; 
   NDSolveValue[{-Laplacian[u[x, y], {x, y}] == f[x, y] 
      + NeumannValue[0, True] // Evaluate,DirichletCondition[u[x, y] == 0, True]}, 
    u[x, y], {x, y} \[Element] mesh0],{i, Length[idx]}];

Then I can plot the basis

Plot3D[basis[[;; ;; 5]], {x, y} \[Element] mesh0, 
 PlotStyle -> Opacity[0.4], PlotRange -> All, PlotTheme -> "Mesh"]

Of course the main point of using the mesh of the FEM is that it can be non trivial. For instance

Needs["NDSolve`FEM`"];
mesh0 = ToElementMesh[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]], 
MaxCellMeasure -> 0.25, AccuracyGoal -> 1]; mesh0["Wireframe"]

while the same code will work exactly unchanged

pol = Table[mesh0["Coordinates"][[ idx[[i]]]] // ConvexHullMesh, {i,Length[idx]}];  
basis = Table[f[x_, y_] = Boole[ RegionMember[pol[[i]], {x, y}]]; 
   NDSolveValue[{-Laplacian[u[x, y], {x, y}] == f[x, y] + 
        NeumannValue[0, True] // Evaluate,
     DirichletCondition[u[x, y] == 0, True]}, 
    u[x, y], {x, y} \[Element] mesh0],{i, Length[idx]}];

And once again

Plot3D[basis[[;; ;; 5]], {x, y} \[Element] mesh0, 
 PlotStyle -> Opacity[0.4], PlotRange -> All, PlotTheme -> "Mesh"]

Now the inverse problem is quite simple

I find the FEM toolkit extremely useful in this context because building a basis function for non trivial geometry is ... non trivial, while the FEM package takes care of it all here.

This solution does not fully address my original question because the basis are 0-splines. Ideally cubic spline would be good too.

Proof of concept for inversion

Let's see how the basis can be used to fit some model. Let us start with a basis on the mesh

basis0 = Table[Boole[ RegionMember[pol[[i]], {x, y}]], {i,Length[idx]}];

and some add hoc density

source[x_, y_] = basis0[[{150, 170, 125}]].{2, 4, 5};
 ContourPlot[source[x, y], {x, y} \[Element] mesh0, PlotPoints -> 75, 
 ContourShading -> None]

that we will try and recover with the corresponding potential:

sol0 = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 
      source[x, y] + NeumannValue[0, True] // Evaluate,
    DirichletCondition[u[x, y] == 0, True]},  u, {x, y} \[Element] mesh0];
Plot3D[sol0[x, y], {x, y} \[Element] mesh0, PlotStyle -> Opacity[0.4],
  PlotRange -> All, PlotTheme -> "ZMesh", PlotPoints -> 50]

Let us sample this potential on a set of random points

data0 = RandomPoint[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]],500] // Sort;

ListPlot[data0, AspectRatio -> 1]

and build the corresponding data set with the value of the potential on those points

data = Map[{#[[1]], #[[2]], sol0[#[[1]], #[[2]]]} &, data0];

Then the linear model follows from the knowledge of the data, data and the basis basis:

ff = Function[{x, y}, basis // Evaluate];
a = Map[ff @@ # &, Most /@ data];
a//Image//ImageAjust 

(looks a bit like the matrix) and we can fit the data as

fit[x_, y_] = basis.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)];

which is a pretty good fit!

Plot3D[fit[x, y] - sol0[x, y], {x, y} \[Element] mesh0,PlotRange -> All]

Similarly we can invert for the source density

inv[x_, y_] =basis0.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)];
Plot3D[inv[x, y], {x, y} \[Element] mesh0, PlotRange -> All, 
 PlotTheme -> "ZMesh", PlotStyle -> Opacity[0.6]]

Of course this inversion is a bit of an overkill to just get the density from the known potential BUT the framework works for any boundary condition and any sampling and arbitrary PDEs that mathematica can solve using FEM.

Indeed, compared to the analytic B-spline method, no extra work in needed to match the boundary conditions since the Mesh generator and FEM package takes care of that.

It is also worth pointing out that once a is known any data set can be inverted almost instantaneously.

To Do

1. I would be best to be able to define cubic splines as well on the mesh (following e.g. this).
2. One needs to write regularisation matrices on the mesh as well, in order to be able to invert ill-conditioned problems (but see this).

Answer 3

Thanks to @Henrik Schumacher's great help in extracting linear piecewise elements from FEM, let me provide a 1-spline solution appropriate for April's fool day.

2D case

Let us start with a fish implicit equation.

reg = ImplicitRegion[(2 x^2 + y^2)^2 - 2 Sqrt[1] x ( 2 x^2 - 3 y^2) + 2 (y^2 - x^2)<= 0, {x, y}]

and discretise it

R = ToElementMesh[R0=DiscretizeRegion[reg], MaxCellMeasure -> 0.015, 
"MeshOrder" -> 1, MeshQualityGoal ->1]; R0

pts = R["Coordinates"]; n = Length[pts];
vd = NDSolve`VariableData[
     {"DependentVariables","Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {R}];
cdata = InitializePDECoefficients[vd, sd,"DiffusionCoefficients" ->
      {{-IdentityMatrix[1]}}, "MassCoefficients" -> {{1}}];
mdata = InitializePDEMethodData[vd, sd];

Discretisation yields

dpde = DiscretizePDE[cdata, mdata, sd];
stiffness = dpde["StiffnessMatrix"];
mass = dpde["MassMatrix"];

To see how it works, let us excite one basis element close to coordinate (0.4,0.1)

i = Nearest[pts -> "Index", {0.4, 0.1}][[2]];
hatfun = ConstantArray[0., n];hatfun[[i]] = 1.;

This is how to interpolate it.

hatfuninterpolated = ElementMeshInterpolation[{R}, hatfun];
plot1 = Plot3D[hatfuninterpolated[x, y], {x, y} \[Element] R, 
  NormalsFunction -> None, PlotPoints -> 50, PlotTheme -> "Business", 
  BoxRatios -> {2, 1, 1}]

In order to compute the corresponding potential let us extract the systemmatrix

bndplist = 
  Sort@DeleteDuplicates[Flatten[R["BoundaryElements"][[All, 1]]]];
intplist = Complement[Range[n], bndplist];

This is what DeployBoundaryConditions does to the stiffness matrix

systemmatrix = stiffness;
systemmatrix[[bndplist]] = 
  IdentityMatrix[n, SparseArray, 
    WorkingPrecision -> MachinePrecision][[bndplist]];

Factorizing the system matrix:

S = LinearSolve[systemmatrix, Method -> "Pardiso"];
load = mass.hatfun;

Solving the actual equation yields the potential for this basis element.

solution = S[load];
solutioninterpolated = ElementMeshInterpolation[{R}, solution];
plot1 = Plot3D[solutioninterpolated[x, y] // Evaluate, 
 {x, y} \[Element] R,NormalsFunction -> None, PlotRange -> All, 
  ColorFunction -> 
   Function[{x, y, z}, RGBColor[1 - z/2, 1 - z, 1/2 - z]], 
  PlotTheme -> "Business", BoxRatios -> {2, 1, 1}]

Let us now define a basis function

basis0 = Table[
   hatfun = ConstantArray[0., n];
   hatfun[[i]] = 1;
   ElementMeshInterpolation[{R}, hatfun],
   {i, 1, n}];

and compute its image

basis = Table[hatfun = ConstantArray[0., n];
   hatfun[[i]] = 1; load = mass.hatfun;solution = S[load];
  ElementMeshInterpolation[{R}, solution],
   {i, 1, n}];

Let us now pick a set of points on our fish

data0 = RandomPoint[R0, 1500] // Sort;
ListPlot[data0]

and define a 'measured potential' from an (ad hoc random) set of 50 basis elements

hatfun0 = ConstantArray[0., n];
hatfun0[[RandomChoice[Range[n], 50]]] = 1;
load = mass.hatfun0;
solution = S[load];
sol0 = ElementMeshInterpolation[{R}, solution];
data = Map[{#[[2]], #[[1]], sol0[#[[2]], #[[1]]]} &, data0];

The linear model relating the basis to the data reads

ff = Function[{x, y}, Map[#[x, y] &, basis] // Evaluate];
a = Map[ff @@ # &, Most /@ data];

Clear[fit];
fit[x_, y_] := Module[{vec = Map[#[x, y] &, basis]},
   vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];

Let us plot the fit:

Plot3D[fit[x, y] // Evaluate, {x, y} \[Element] R, 
 NormalsFunction -> None, PlotRange -> All, 
 ColorFunction -> 
  Function[{x, y, z}, RGBColor[1 - z/2, 1 - z, 1/2 - z]], 
 PlotTheme -> "Business", BoxRatios -> {2, 1, 1}]

We can now also invert it:

Clear[inv];
inv[x_, y_] := Module[{vec = Map[#[x, y] &, basis0]},
   vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];
Plot3D[inv[x, y] // Evaluate, {x, y} \[Element] R, 
 NormalsFunction -> None, 
 ColorFunction -> Function[{x, y, z}, 
 RGBColor[1 - z/2, 1 - z, 1/2 - z]], 
 PlotTheme -> "Business", PlotPoints -> 50, BoxRatios -> {2, 1, 1}, 
 PlotRange -> {0, 2}]

It compares well with the input model:

hatfuninterpolated = ElementMeshInterpolation[{R}, hatfun0];
plot1 = Plot3D[hatfuninterpolated[x, y], {x, y} \[Element] R, 
  NormalsFunction -> None, PlotPoints -> 50, PlotTheme -> "Business", 
  BoxRatios -> {2, 1, 1},
  PlotRange -> {0, 2}]

Caveat: this is most likely not as efficient as it should be (see Henrik's comments). I could imagine e.g. that the way the basis function is defined is probably in part redundant w.r.t. to what is available within the FEM toolbox.

But it does illustrate that we can invert a given PDE with arbitrary sampling and ad hoc boundary condition on a set of linear piecewise basis function, which is differentiable, which is pretty cool IMHO. This question/answer provides means of regularising the inversion should this be needed (i.e. if a is poorly conditioned, with very small eigenvalues).

3D case

Let us give in one block the 3D code on a unit ball:

R = ToElementMesh[R0 = Ball[], MaxCellMeasure -> 0.125/16, 
AccuracyGoal -> 1, "MeshOrder" -> 1];pts = R["Coordinates"];n = Length[pts];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y, z}}];
sd = NDSolve`SolutionData[{"Space"} -> {R}];
cdata = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-IdentityMatrix[3]}}, 
   "MassCoefficients" -> {{1}}];
mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd];
stiffness = dpde["StiffnessMatrix"];
mass = dpde["MassMatrix"];
bndplist = Sort@DeleteDuplicates[Flatten[R["BoundaryElements"][[All, 1]]]];
intplist = Complement[Range[n], bndplist]; systemmatrix = stiffness;
systemmatrix[[bndplist]] = 
  IdentityMatrix[n, SparseArray, 
    WorkingPrecision -> MachinePrecision][[bndplist]];
S = LinearSolve[systemmatrix, Method -> "Pardiso"];
   basis0 = Table[
   hatfun = ConstantArray[0., n];
   hatfun[[i]] = 1;
   ElementMeshInterpolation[{R}, hatfun],
   {i, 1, n}];
   basis = Table[
   hatfun = ConstantArray[0., n];
   hatfun[[i]] = 1; load = mass.hatfun;
   solution = S[load];
   solutioninterpolated = ElementMeshInterpolation[{R}, solution];
   solutioninterpolated,
   {i, 1, n}];

data0 = RandomPoint[R0, 500] // Sort;    
hatfun0 = ConstantArray[0., n];
hatfun0[[RandomChoice[Range[n], 50]]] = 1;
load = mass.hatfun0; solution = S[load];
sol0 = ElementMeshInterpolation[{R}, solution];

data = Map[{#[[1]],#[[2]],#[[3]],sol0[#[[1]], #[[2]],#[[3]]]} &, data0];
ff = Function[{x, y, z}, Map[#[x, y, z] &, basis] // Evaluate];
a = Map[ff @@ # &, Most /@ data];   
Clear[fit];
fit[x_, y_, z_] := Module[{vec = Map[#[x, y, z] &, basis]},
   vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];  
Clear[inv];
inv[x_, y_, z_] := Module[{vec = Map[#[x, y, z] &, basis0]},
   vec.LinearSolve[Transpose[a].a, Transpose[a].(Last /@ data)]];

As a check let us look at the cross section through the mid-plane of the inverted density and the input density resp.

Plot3D[inv[x, y, 0] // Evaluate, {x, y} \[Element] Disk[], 
 NormalsFunction -> None, ColorFunction -> 
  Function[{x, y, z}, RGBColor[1 - z/2, 1 - z, 1/2 - z]], 
 PlotTheme -> "Business", PlotPoints -> 50, BoxRatios -> {1, 1, 1}, 
 PlotRange -> {0, 2}]

hatfuninterpolated = ElementMeshInterpolation[{R}, hatfun0];
plot1 = Plot3D[hatfuninterpolated[x, y, 0], {x, y} \[Element] Disk[], 
  NormalsFunction -> None, PlotPoints -> 50, PlotTheme -> "Business", 
  BoxRatios -> {1, 1, 1},PlotRange -> {0, 2}]

It seems to work fine!