# Piecewise Continuous Linear Basis on a Triangular Mesh

Question

Given a FEM mesh, I would like to define a set of basis functions anchored to the mesh, so that any piecewise linear continuous function on the mesh can be expanded over that set.

Such a basis is possibly called ElementShapeFunction in the FEM framework?

Attempt

Let us consider

mesh0 = ToElementMesh[RegionUnion[Disk[], Rectangle[{0, 0}, {2, 2}]],
MaxCellMeasure -> 0.125, AccuracyGoal -> 1,
MeshQualityGoal -> "Maximal",
"MeshOrder" -> 1, "MeshElementType" -> TriangleElement]

Given this mesh, I should be able to construct a set linear piecewise generalisation of indicator function (pyramide-like) which should look like this:

where the top dark red vertex is above the middle red vertex.

I understand that the FEM tools allow me to identify the triangles

me = mesh0["MeshElements"][[1, 1]]; nn = Length[me];
Table[{Hue[i/nn], Polygon@mesh0["Coordinates"][[me[[i]]]]}, {i, nn}] // Graphics

Thanks to @user21, I also gathered how to extract the linear piecewise function on a regular (unit) triangle

Table[Table[{r, s,
ElementShapeFunction[TriangleElement, 1][r, s][[i]]},
{s, 0, 1, 1/100}, {r, 0, s, 1/100}] // Flatten[#, 1] & //
ListContourPlot[#, Axes -> True] &,
{i, 3}]

I also vaguely recall that the transformation to the regular triangle involve the inverse of $$\left( \begin{array}{ccc} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \\ \end{array} \right)$$ where the $$(x_i,y_i)$$ are the coordinates of the 3 vertices of the triangles.

So in principle I am all set to define my basis function.

But my target to have a basis which is consistent with the way things are done within the FEM package. This I cannot easily do with some help from people who know its internals. I truly believe this basis could be useful to many projects, beyond my own, especially if it is consistent with the logic of the package.

So I am after a function which would take mesh as an argument and return a list of $$\cal N_i$$ functions such as plotted above, so that any linear function on the mesh can be written unambiguously as a sum over these.

The format of these functions should ideally generalise that of the existing BSplineFunction for triangular meshes. Internally they should correspond to a piecewise description of the linear interpolation over the relevant triangles.

1. One complication I can see arises from e.g. the top left blue triangle which is not part of a polyhedra.
2. Ideally one wants to have special basis elements on the edges which can be non zero, or with a given slope.

The latter requirement is possibly achieved by counting the outer vertex as a double (or triple) knot, as is done for BSplines?

Eventually, this basis could be replacing the constant piecewise function, or for a regular mesh those presented in this answer .

• The point to efficiency in the finite element method is to avoid explicit representations of the basis functions, but to work only with the coordinates with respect to such a basis. If you generate your mesh with ToElementMesh with MeshOrder->1 then mass and stiffness matrix are already assembled with respect to this basis. If you have $n$ vertices, then the "hat function" of vertex $i$ ($=1$ at that vertex, $=0$ at all other vertices and piecewise linear in between) is represented by the i-th basis vector of the standard basis on $\mathbb{R}^n$. – Henrik Schumacher Mar 31 at 15:55
• To put that hat function on "the right hand side" of the equation (and that's what I guess is what you are about to do), it typically suffices to multiply the corresponding vector my the mass matrix (and to adjust the boundary degrees of freedom to match the boundary conditions). At least, that's the case for typical elliptic equations. – Henrik Schumacher Mar 31 at 15:56
• @HenrikSchumacher you know better than me, BUT it seems to me that mathematica has both LinearInterpolation and BSplineBasis which serve different purpose. And I am trying to implement 1-Spline on a triangular mesh because you showed me how to write Laplacian Penalty function on such basis :-) – chris Mar 31 at 16:00
• @HenrikSchumacher re':it typically suffices to multiply' please please please could you give an example while answering this question mathematica.stackexchange.com/questions/216747/…? I apologise to be so clueless about FEM. – chris Mar 31 at 16:01
• Even if it is inefficient it could be educationally useful to be able to access such basis? In the spirit of a top down approach to learning FEM (as discussed by @user21 at some earlier stage on this site), or in the spirit of this elegant answer mathematica.stackexchange.com/a/110210/1089 – chris Mar 31 at 16:24

## 1 Answer

Okay, here a small example for an elliptic boundary value problem.

(I am reusing some old code over and over again, so user21 will certainly remind me again that a couple of things can be simplified... ;) )

Let's start with my favorite region.

Needs["NDSolveFEM"]
(*Initialization of Finite Element Method*)
R = ToElementMesh[
BoundaryMeshRegion[
Map[t \[Function] (2 + Cos[5 t])/3 {Cos[t], Sin[t]},
Most@Subdivide[0., 2. Pi, 2000]],
Line[Partition[Range[2000], 2, 1, 1]]
],
MaxCellMeasure -> 0.001,
"MeshOrder" -> 1
];
pts = R["Coordinates"];
n = Length[pts];
vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {R}];
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-IdentityMatrix[2]}},
"MassCoefficients" -> {{1}}
];
mdata = InitializePDEMethodData[vd, sd];

(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd];
stiffness = dpde["StiffnessMatrix"];
mass = dpde["MassMatrix"];

This supplies us with a stiffness matrix stiffness and a mass matrix mass, both assembled with repect to a basis of piecewise linear hat functions. We will need them only

Let's choose a vertex somewhere in the middle and represent its hat function as a vector with respect to the basis function. (In the following, I put my comments into the codes so that it is easier to copy.)

i = Nearest[pts -> "Index", {0., 0.1}][[1]];
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]; // AbsoluteTiming // First

(*But the interpolation is actually not needed because the graph of the function can be plotted like this:*)
scale = 2/3;
plot2 = Graphics3D[{
GraphicsComplex[Join[pts, scale Partition[hatfun, 1], 2],
Polygon[R["MeshElements"][[1, 1]]]]
}]; // AbsoluteTiming // First

GraphicsRow[{plot1, plot2}, ImageSize -> Large]

0.251001

0.000127

Notice the difference in the timings. That's basically the reason why I say you should avoid interpolation function as much as possible.

Okay, let's go an. We want to see how to use hatfun as the right hand side $$b$$ of the pde $$\Delta u = b$$ in $$\varOmega$$ and $$u|_{\partial \varOmega} = f$$.

(*Finding boundary and interior degrees of freedoms.*)

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"];

This is all that we have to do for the system matrix.

(*This is how the NDSolveFEM builds the load vector (a.k.a. the \
right hand side). *)
load = mass.hatfun;

(*f is a function that specifies the Dirichlet boundary conditions.*)

f = {x, y} \[Function] 0.0001 Sin[25 ArcTan[x, y]];
(*This is what DeployBoundaryConditions does to the load vector*)

load[[bndplist]] = f @@@ pts[[bndplist]];

(*Solving the actual equation.*)
solution = S[load];

(*Plotting via interpolation.*)

solutioninterpolated = ElementMeshInterpolation[{R}, solution];
plot1 = Plot3D[solutioninterpolated[x, y], {x, y} \[Element] R,
NormalsFunction -> None, PlotRange -> All]; //
AbsoluteTiming // First

(*Fast plotting*)
scale = 1200;
plot2 = Graphics3D[{
GraphicsComplex[Join[pts, scale Partition[solution, 1], 2],
Polygon[R["MeshElements"][[1, 1]]]]
}]; // AbsoluteTiming // First

GraphicsRow[{plot1, plot2}, ImageSize -> Large]

0.241259

0.000119