Penalty function on discrete mesh using Laplace-Beltrami operator?

Question

Context

I am interested in extending to the ill-condionned regime the inversion of linear equations arising from inverting differential equations which have been solved via 0-splines over a mesh using the FEM toolkit in connection to this question.

For this purpose I need to compute a discrete Laplacian operator on a given mesh produced by ToElementMesh.

Question

Given some mesh, and a discrete function associating a value at each mesh element, I would like to compute a penalty function corresponding to the integral of the Laplacian square of the function over that mesh.

$$ P(\mathbf{a}) = \int \big|\Delta \phi \big|^2 dx\,, $$

where $\mathbf{a}=({a_i}_{i\le n})$ is a vector of values on the mesh elements and $\phi(\mathbf{x})=\sum_i a_i \phi_i(\mathbf{x})$, with $\phi_i(\mathbf{x})=1$ iff $\mathbf{x}\in \mbox{cell}_i$ and $0$ otherwise.

Of course strictly speaking, as defined, $P$ is formally zero almost everywhere since the functions $\phi_i$ are constant.

What I am after is a Sparse matrix, $\cal D$, so that

$$P(\mathbf{a}) = \mathbf{a}^T\cdot \cal D \cdot \mathbf{a}. $$

I am fairly sure some element of the answer is available in the amazing answer involving the Laplace-Beltrami operator.

It would be best if the answer would work with meshes in dimension 2 and 3.

Attempt

I have implemented a test case.

 mesh0 = ToElementMesh[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]], 
  MaxCellMeasure -> 0.125, AccuracyGoal -> 2]
 mesh0["Wireframe"]

From the mesh I can find their centroid

idx = mesh0["MeshElements"][[1, 1]];
tt = Table[mesh0["Coordinates"][[ idx[[i]]]], {i, Length[idx]}];
center = Map[1/Length[#] Plus @@ # &, tt, {1}];
ListPlot[center, AspectRatio -> 1]

I can then compute the matrix of distances between centroids of the mesh elements

 dist = DistanceMatrix[center]; 

If my mesh was regular I could use

s = SparseArray[{{i_, i_} -> -1, {i_, j_} /; i - j == 1 -> 
      2, {i_, j_} /; i - j == 2 -> -1}, {17, 15}] // Transpose;
s1 = ArrayFlatten[TensorProduct[s, s]];
pen = Transpose[s1].s1; pen // MatrixPlot

So an alternative is to compute difference of values at 3 centres, $2x_i -x_{i-1}-x_{i+1}$ and divide by the distance square between those centres as a discrete proxy for the Laplacian.

dif = SparseArray[{{nn, nn} -> 1,
    {1, 1} -> 1, {i_, i_} -> 
     2, {i_, j_} /; i - j == 1 -> -1, {i_, j_} /; 
      i - j == -1 -> -1}, {nn, nn}];
idist = Inverse@DistanceMatrix[center] // SparseArray;
idist = Transpose[idist]. idiot;
pen = Transpose[idist.dif].(idist.dif); pen // MatrixPlot

This operator has the good taste to nul a constant vector, but it is expansive to compute. May be a workaround with Nearest is in order to make idist sparse?

Answer 1

If you insist on using my code from the other post, this can be done as follows; note that you have to settle with first order FEM in that case:

Needs["NDSolve`FEM`"];

mesh0 = ToElementMesh[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]],
  MaxCellMeasure -> 0.125, AccuracyGoal -> 2, "MeshOrder" -> 1]
pts = mesh0["Coordinates"];
faces = mesh0["MeshElements"][[1, 1]];
pat = Partition[Flatten[getSurfaceLaplacianCombinatorics[faces]], 2];
flist = Flatten[faces];
laplacian = SurfaceLaplaceBeltrami[pts, flist, pat];
mass = SurfaceMassMatrix[pts, flist, pat];
lumpedmass = Total[mass];
invlumpedmass = 
  SparseArray[
    Partition[Union @@ mesh0["BoundaryElements"][[1, 1]], 1] -> 0., 
    Length[mass], 1.]/Total[mass];
bilaplacian = laplacian.(invlumpedmass laplacian);

Then the matrix bilaplacian is a somewhat crude approximation to what you are looking for. It is crude because it uses nonconforming (first order) finite elements and because it uses mass lumping. But it should be quite fast because of this (inverting the mass matrix mass would lead to a dense matrix). You only require it as a regularizer, so this should work out well.

In general, you can use the stiffness matrix laplacian and the mass matrix mass provided by any other FEM tool (e.g., you can obtain them from the low level FEM tools in Mathematica, too). The only other ingredient would be a diagonal matrix A with ones on the diagonal for interior degrees of freedom and zeroes for the boundary degrees of freedom. Then the matrix that you seek should be

bilaplacian = laplacian.A.Inverse[mass].A.laplacian

Typically, Inverse[mass] is a dense matrix, so one should avoid inverting mass if possible. With first order FEM, one can employ mass lumping (as I did above). From what I heard, mass lumping does not work well for higher order FEM (but I could be wrong). Hence I would suggest Mathematica first order low level FEM tools for the 3D case. For the 2D case with a planar mesh, it is up to you which one you want to use. I do not know whether Mathematica supports surface FEM in version 12.1; it does not in version 12. So if you want to use that for surfaces, you are doomed to use my code, I guess. ;)

Answer 2

Following @HenrikSchumacher's and @user21's advice, I have extracted from the FEM tutorial a computation of stiffness matrix on the mesh as follows

Needs["NDSolve`FEM`"];
mesh = ToElementMesh[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]], 
       MaxCellMeasure -> 0.125, AccuracyGoal -> 1, "MeshOrder" -> 1];
 nr = ToNumericalRegion[mesh];
coefficients={"DiffusionCoefficients"->{{IdentityMatrix[2]}},"LoadCoefficients"->{{1}}};
vd = NDSolve`VariableData[{"DependentVariables" -> {u},"Space" -> {x, y}}];
sd = NDSolve`SolutionData[{"Space" -> nr}];
initCoeffs = InitializePDECoefficients[vd, sd, coefficients];
methodData = InitializePDEMethodData[vd, sd];
finiteElements=DiscretizePDE[initCoeffs,methodData, sd,"SaveFiniteElements" -> True];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];

Then

    stiffness // MatrixPlot

is (hopefully!) a matrix which applies a Laplacian to the coefficients of the (piecewise linear) 1-spline evaluated on the vertices of the mesh.

Nicely, the method works for 3D meshes as well

Needs["NDSolve`FEM`"];
mesh = ToElementMesh[Ball[],MaxCellMeasure->0.125/8,AccuracyGoal->1, "MeshOrder" -> 1];
Show[{mesh["Wireframe"],  mesh["Coordinates"] // 
   ListPointPlot3D[#, PlotStyle -> AbsolutePointSize[10],
     ColorFunction -> Function[{x, y, z}, RGBColor[x, y, z]]] &}]

nr = ToNumericalRegion[mesh];
coefficients = {"DiffusionCoefficients" -> {{IdentityMatrix[3]}}, 
   "LoadCoefficients" -> {{1}}};
vd = NDSolve`VariableData[{"DependentVariables" -> {u}, 
    "Space" -> {x, y, z}}];
sd = NDSolve`SolutionData[{"Space" -> nr}];
initCoeffs = InitializePDECoefficients[vd, sd, coefficients];
methodData = InitializePDEMethodData[vd, sd];
finiteElements=DiscretizePDE[initCoeffs, methodData, sd,"SaveFiniteElements" -> True];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
stiffness // MatrixPlot