The Variational Finite Element Method for Laplace and Poission Equations

Section 1: Setting up the equations for the variational method.

The basic idea of the variational method is that the solution to 

Δu = r

is also the minimizer of the functional

D(u,u) = ∫ (1/2) ((‖ ∇u ‖)^2) + ru dVol

If we triangulate our domain, then we can define a piecewise linear function on each triangle. The Laplacian is a little hard to define on such a function (where is the second derivative defined?), but the functional D(u,u) makes perfect sense. Since integration is linear we can then write down the total integral as a sum of contributions from each triangle.

Assume that a given triangle is specified by ((Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])), ((Subscript[x, j]),(Subscript[y, j]),(Subscript[z, j])), ((Subscript[x, k]),(Subscript[y, k]),(Subscript[z, k])), and we want to write down the gradient of the linear function through these points. We can rewrite this function as a sum of standard basis functions

(Subscript[L, i]) = plane through ((Subscript[x, i]),(Subscript[y, i]),1), ((Subscript[x, j]),(Subscript[y, j]),0), ((Subscript[x, k]),(Subscript[y, k]),0)
(Subscript[L, j]) = plane through ((Subscript[x, i]),(Subscript[y, i]),0), ((Subscript[x, j]),(Subscript[y, j]),1), ((Subscript[x, k]),(Subscript[y, k]),0)
(Subscript[L, k]) = plane through ((Subscript[x, i]),(Subscript[y, i]),0), ((Subscript[x, j]),(Subscript[y, j]),0), ((Subscript[x, k]),(Subscript[y, k]),1)

We now compute the gradient of (Subscript[L, i]).

Li = Collect[Simplify[InterpolatingPolynomial[{{{xi,yi},1},{{xj,yj},0},{{xk,yk},0}},{X,Y}]],{X,Y}]
GradLi = Outer[D,{Li},{X,Y}][[1]]

(Subscript[L, j]) and (Subscript[L, k]) are similar:

Lj = Collect[Simplify[InterpolatingPolynomial[{{{xi,yi},0},{{xj,yj},1},{{xk,yk},0}},{X,Y}]],{X,Y}]
GradLj = Outer[D,{Lj},{X,Y}][[1]]

Lk = Collect[Simplify[InterpolatingPolynomial[{{{xi,yi},0},{{xj,yj},0},{{xk,yk},1}},{X,Y}]],{X,Y}]
GradLk = Outer[D,{Lk},{X,Y}][[1]]

SquaredNorm = Simplify[Dot[zi * GradLi + zj * GradLj + zk*GradLk,zi * GradLi + zj * GradLj + zk*GradLk]]

Now the point here is that these are linear functions. So their gradients are constant across the triangle. Which means that integrating the squared norm of the gradient over the triangle is the same as just multiplying the SquaredNorm function by the area of the triangle.

Area = (1/2)*Cross[{xj,yj,0} - {xi,yi,0},{xk,yk,0} - {xi,yi,0}][[3]]

We end up with the following norm function, which collects the entire value of 

 ∫ (1/2) ((‖ ∇u ‖)^2) dVol
 
 across this triangle:

NormFun = (1/2) * Collect[SquaredNorm * Area,{zi,zj,zk}]

This is a quadratic polynomial in the function variables {zi,zj,zk} which depends on the mesh variables {xi,xj,xk}, {yi,yj,yk} as well, as we can see by using the CoefficientRules function to write down the whole thing as a multinomial in zi, zj, zk. 

Here the first triples are powers of the three variables, so {1,1,0} means zi^1 zj^2 zk^0, while the thing to the right of the rule is the coefficient, which you can see is purely a function of the mesh geometry and doesn’t actually involve the zi, zj, zk variables at all!

Simplify[CoefficientRules[NormFun,{zi,zj,zk}]]

We now switch to considering the term 

∫rudVol

in our expression for D(u,u). In this term, we’ll use a linear approximation to the r function given by values at the mesh coordinates. We’ll actually have to integrate this linear function over the triangle, because it’s not a constant function (unlike the ((‖ ∇u ‖)^2) term).

Ri = Simplify[InterpolatingPolynomial[{{{xi,yi},ri},{{xj,yj},rj},{{xk,yk},rk}},{X,Y}]]

RZFun = Ri * (zi * Li + zj*Lj + zk*Lk)

We need to change coordinates to integrate over the triangle. The idea is that if we make the following change of coordinates, the triangle is easy to express in (u,v) coordinates:

RZFunUV = RZFun /. {X -> u(xj - xi) + v(xk - xi), Y -> u (yj - yi) + v(yk - yi)}

And then we can integrate:

RZIntegral = Integrate[RZFunUV,{u,0,1},{v,0,1-u}]

We have to multiply by the Jacobian of the transformation from (u,v) to (x,y) coordinates in order to convet back to an integral in the original dx dy area form.

RZJacobian = Det[Outer[D,{u(xj - xi) + v(xk - xi),u (yj - yi) + v(yk - yi)},{u,v}]]

And we end up, finally, with the following (fairly grim) term in the {ri, rj, rk} numbers and {zi,zj,zk} variables. Note that this is linear in {zi, zj, zk}.

RZTerm = Collect[Simplify[RZJacobian * RZIntegral],{zi,zj,zk}]

Now to find a minimizer of this problem, we can use the least-squares picture from class:



This tells us that our functional (which is after all just distance from the origin in some moral sense) is a quadratic function which should have a single critical point, which is the global minimum of the function. This means that we can solve the (quadratic) minimization problem by computing the gradient vector of the functional 

D(u,u) = ∫ (1/2) ((‖ ∇u ‖)^2) + ru dVol

with respect to the zi variables (which, remember, are the only actual variables in this problem-- everything else just specifies our choice of mesh), setting 

grad_z D(u,u) = 0

and viewing that as a large, sparse linear system of equations in the zi which we can solve numerically. These are called the Galerkin equations.

We’re first going to compute the contribution to the partial of D(u,u) with respect to zi, zj, and zk coming from this triangle. Notice that this is not an entire Galerkin equation. To get the whole thing, we’ll have to sum up these contributions for every triangle where vertex i appears.

So let’s begin by differentiating the equations for NormFun on our triangle with respect to zi, zj, and zk. This gives the first part of our equations (this is the linear part in zi, zj, and zk).

NormGalerkinEquations = Outer[D,{NormFun},{zi,zj,zk}][[1]]

As you might expect, it’s a little easier to see that these are actually linear in the zi, zj, and zk if we go ahead and rewrite the coefficients explicitly using simplify and CoefficientRules.

Simplify[CoefficientRules[#,{zi,zj,zk}]] & /@ NormGalerkinEquations

Now what about that RZ term? 

∫ ru dVol

Well, we can differentiate it with respect to the zi, zj, and zk variables, too. We get a collection of three constants:

RZGalerkinEquations =Simplify /@ Outer[D,{RZTerm},{zi,zj,zk}][[1]]

Again, checking CoefficientRules verifies that these don’t contain zi, zj, or zk terms:

Simplify[CoefficientRules[RZGalerkinEquations,{zi,zj,zk}]]

We’ll see shortly that things change when we have one or more boundary nodes involved. So for the moment, we simply collect the entire (quadratic) objective function into a single term:

ObjectiveFunction = NormFun + RZTerm

Section 2: Mesh generation

Now, of course, we need to triangulate our domain! This is really not an easy problem, and in fact mesh generation is a current research problem in applied math. For the moment, we’re going to punt and load some pregenerated meshes. A good 2d option is still Jonathan Shewchuk’s triangle (old, but reliable). And Mathematica includes a 3d mesh generator called TetGen.

LakeNodes =  Import[NotebookDirectory[] <> "lake_nodes.txt","table"];
LakeTriangles = Import[NotebookDirectory[] <> "lake_elements.txt","table"];

MakeGCAndBoundaryList[Nodes_,Triangles_] := Module[{Edges,EdgeComplex,TriComplex,BoundaryEdges,BoundaryComplex,BoundaryNodes,TrimmedTriangles},

Edges= Sort /@ Flatten[Subsets[#,{2}] & /@ Triangles,1];
BoundaryEdges = #[[1]] & /@ Select[Tally[Edges],#[[2]] == 1&];
BoundaryNodes = Union[Flatten[BoundaryEdges]];
TrimmedTriangles = Select[Triangles,Length[# ⋂ BoundaryNodes] ≠ 3&];

TriComplex = {EdgeForm[Opacity[0.2]],FaceForm[ColorData["AlpineColors"][RandomVariate[UniformDistribution[{0.8,1}]]]],#}& /@(Normal[GraphicsComplex[Nodes, Polygon /@ TrimmedTriangles]]);

Edges= Sort /@ Flatten[Subsets[#,{2}] & /@ TrimmedTriangles,1];
BoundaryEdges = #[[1]] & /@ Select[Tally[Edges],#[[2]] == 1&];

BoundaryComplex = {ColorData["AlpineColors"][0.5],GraphicsComplex[Nodes,Line /@ BoundaryEdges]};

{Graphics[{TriComplex,BoundaryComplex}],BoundaryNodes,Sort[Complement[Range[Length[Nodes]],BoundaryNodes]],TrimmedTriangles}
];

There is one clever piece to this code which is worth a little explanation:

Tally counts the number of times each element occurs in a list. If an edge occurs twice (in two different triangles) it is an interior edge-- otherwise, it occurs once and is a boundary edge.

Further, we take the liberty of deleting any triangle where all of the nodes are boundary nodes. This triangle can’t affect the solution anyway because the function is entirely fixed over this triangle, and it may confuse things later on.

Note that we have to recompute the boundaryedges afterwards to delete any boundary nodes which only appeared on a triangle fully contained in the boundary.

{LakeImage,LakeBoundary, LakeVars,LakeTrimmedTriangles} = MakeGCAndBoundaryList[LakeNodes,LakeTriangles];

LakeImage

LakeBoundary

LakeVars

Length[LakeVars]

Length[LakeBoundary]

Section 3. Assembling the Galerkin Matrix

We now assume that we’re going to have Dirichlet boundary conditions with u = 0 on the boundary and a mesh given by a list of pairs of coordinates and a list of triples of coordinate indices forming triangles. Further, we’re going to get a list of boundary vertices. Our job is to separate the actual variables from the boundary variables and generate a SparseArray which we can use for the linear algebra. We’re going to have work out the r values from an rFunction, which we expect to be a scalar-valued function R[x_,y_]. 

We’re going to write down the equations corresponding to 

∇(objective function) = 0

keeping in mind that the objective function is a sum of triangle functions and that the variables are given by whichever of the zi coordinates are NOT on the boundary.

MakeTriangleMatrix[Vars_,Boundary_,tri_,Nodes_,R_,ObjectiveFunction_] := Module[{fun,trules,boundaryrules,zvars,zZeroRules,b,A},

" The first thing to do is for each triangle, substitute in everything we know (the xi, yi, xj, yj, xk, yk coordinates, the ri, rj, rk values, and zeros for any of zi, zj, or zk which are on the boundary).";

trules = {xi->Nodes[[tri[[1]]]][[1]],yi -> Nodes[[tri[[1]]]][[2]],xj->Nodes[[tri[[2]]]][[1]],yj -> Nodes[[tri[[2]]]][[2]],
xk -> Nodes[[tri[[3]]]][[1]],yk->Nodes[[tri[[3]]]][[2]],
ri -> R[xi,yi], rj -> R[xj,yj], rk -> R[xk,yk]};
boundaryrules = Pick[{zi->0,zj->0,zk->0},MemberQ[Boundary,#]& /@ tri];
fun = ObjectiveFunction //. Join[trules,boundaryrules];
(*Print[fun];*)

" The function 'fun' is now quadratic in whichever of the zi, zj, zk are not on the boundary (and so are actual variables). We now make a list of pairs of variable numbers (positions in the Vars array) and variable names (zi, zj, zk) for surviving variables";

(*Print[Vars];
Print[tri];
Print[MemberQ[Vars,#] & /@ tri];*)

zvars = 
{#[[1]],Position[Vars,#[[2]]][[1,1]],#[[2]]}& /@Pick[Transpose[{{zi,zj,zk},tri}],MemberQ[Vars,#] & /@ tri];

zZeroRules = #[[1]] -> 0 & /@ zvars; 
(*Print[zvars];*)

" The contributions to the b vector are easy: they are the constant terms in the derivatives of fun with respect to whichever subset of {zi, zj, zk} appear in the zvars list (because they are variables). 

	 To get the constant term, we'll use the zZeroRules, which set each of the surviving members of {zi,zj,zk} to zero.

	 Notice that when we specify the SparseArray for b, we have to give it explicit dimensions (Length[Vars]) and tell Mathematica that we're assuming that entries we don't give here have the value zero. ";

b = SparseArray[#[[2]]->(D[fun,#[[1]]] /. zZeroRules) & /@ zvars,Length[Vars],0];

(*Print[b[[#[[2]]]] & /@ zvars];*)

" The contributions to the A matrix are only a little harder. Here, we need to make pairs of zvars: the {i,j} entry in the A array is the coefficient of variable j in the partial of fun with respect to variable i. 

     Again, we have to specify that the matrix A is Length[Vars] x Length[Vars],
     and that missing entries are taken to be zero in order to be able to sum the SparseArray terms later.";

(*Print[Tuples[{zvars,zvars}]];*)

A = SparseArray[{#[[1,2]],#[[2,2]]} -> Coefficient[D[fun,#[[1,1]]],#[[2,1]]] & /@Tuples[{zvars,zvars}],{Length[Vars],Length[Vars]},0];

{A,b}

];

Since we can sum sparse matrices (and vectors) in Mathematica and get sparse matrices and vectors out, we now just total the contributions over the triangles.

MakeDirichletProblem[Nodes_,Triangles_,Boundary_,Vars_,ObjectiveFunction_,R_] := 
Total /@ (Transpose[MakeTriangleMatrix[Vars,Boundary,#,Nodes,R,ObjectiveFunction]  & /@ Triangles]);

We’re now going to write down the problem for the constant function r = -1.

MyR[x_,y_] := -1.0;

{A,b} = MakeDirichletProblem[LakeNodes,LakeTrimmedTriangles,LakeBoundary,LakeVars,ObjectiveFunction,MyR]

We get a 352 x 352 matrix, which we can plot with MatrixPlot to appreciate the sparsity pattern of the matrix. As you can see, the matrix is kind of loosely banded, and clearly has an interesting structure, but doesn’t have anything really obvious for us to grab onto.

MatrixPlot[A]

Mathematica is perfectly happy to solve a sparse 352x352 linear system almost instantly:

sol = LinearSolve[A,b];

Section 4. Converting the solution vector back into functional data.

We now convert the solution vector back into function data by unwrapping our assignment of vars to generate a 3d graphicscomplex which contains the z-values.

MakeDirichletSolutionGraphic[Nodes_,Triangles_,Boundary_,Vars_,Soln_] := Module[{SpaceNodes,Edges,BoundaryEdges,BoundaryNodes,TrimmedTriangles,TriComplex,BoundaryComplex},

SpaceNodes = MapIndexed[{#1[[1]],#1[[2]],If[MemberQ[Vars,First[#2]],0.01* Soln[[Position[Vars,First[#2]][[1,1]]]],0]} &,Nodes];

(*Print[SpaceNodes];*)

Edges= Sort /@ Flatten[Subsets[#,{2}] & /@ Triangles,1];
BoundaryEdges = #[[1]] & /@ Select[Tally[Edges],#[[2]] == 1&];
BoundaryNodes = Union[Flatten[BoundaryEdges]];
TrimmedTriangles = Select[Triangles,Length[# ⋂ BoundaryNodes] ≠ 3&];

TriComplex = {EdgeForm[Opacity[0.2]],FaceForm[ColorData["AlpineColors"][RandomVariate[UniformDistribution[{0.8,1}]]]],#}& /@(Normal[GraphicsComplex[SpaceNodes, Polygon /@ TrimmedTriangles]]);

Edges= Sort /@ Flatten[Subsets[#,{2}] & /@ TrimmedTriangles,1];
BoundaryEdges = #[[1]] & /@ Select[Tally[Edges],#[[2]] == 1&];

BoundaryComplex = {ColorData["AlpineColors"][0.5],GraphicsComplex[SpaceNodes,Line /@ BoundaryEdges]};

Graphics3D[{TriComplex,BoundaryComplex}]
];

MakeDirichletSolutionGraphic[LakeNodes,LakeTrimmedTriangles,LakeBoundary,LakeVars,sol]

Section 5. Refining the mesh.

Of course, to check this, it would be nice to add different boundary conditions, or change the mesh by replacing each triangle with four inscribed triangles and so forth. I don’t think we’re going to do this for this version of the class, but maybe next time...