Illustrating the Finite Element Method By Solving Laplace' s Equation In Two Dimensions For A Non-Rectangular Region

Frank Kampas (fkampas@msn.com)

Introduction

The finite element method is the standard approach for solving partial differential equations in geometries in which the boundaries are not simple variable ranges.  This article illustrates the basic approach by solving Laplace's equation in two dimensions ( (∂_(x,x))u +(∂_(y,y))u = 0  ) for a region consisting of two unequal rectangles joined together, for a Dirichlet boundary condition (the function is specified on the boundary).

The method used minimizes the square of the gradient of the function u, ∇u · ∇u = ((((∂_x) u))^2)+ ((((∂_y) u))^2) in the region.  The equivalence of this minimization to the 2-d Laplace equation can be seen using the package VariationalMethods.

Needs["VariationalMethods`"]

EulerEquations[((((∂_x) u[x,y]))^2) + ((((∂_y) u[x,y]))^2),u[x,y],{x,y}]

The minimization is performed by dividing the region into triangles and assuming that the solution is a linear function of x and y inside each triangle.  The minimization variables are the values of the function at the nodes of the mesh.

In actual practice, a more complicated function would be used to for the solution and the process would be carried out iteratively, with the elements subdivided repeatedly in order to lower the error to below some threshold.  Despite the simplifications used, the result obtained is in good agreement with a more sophisticated calculation except near the boundary between two closely spaced right-angle bends.

Mathematical Formalism

A linear function of x and y is used inside each triangle, as stated in the introduction.

f[{x_,y_}]=a*x+b*y+c;

The square of the gradient is simply given by (a^2)+(b^2).

((((∂_x) f[{x,y}]))^2) + ((((∂_y) f[{x,y}]))^2)

If the three vertices of the triangle are { x1, y1 }, { x2, y2 } and { x3, y3 } and the values of the function at the vertices are u1, u2 and u3, the coefficients a, b and c can be found as follows:

fsoln=Solve[Thread[f/@ {{x1,y1},{x2,y2},{x3,y3}}=={u1,u2,u3}],{a,b,c}]

The square of the gradient in the triangle, as a function of the coordinates of the vertices and the function values at them,  is

grad2[{{x1_,y1_,u1_},{x2_,y2_,u2_},{x3_,y3_,u3_}}]=((a^2)+(b^2))/. fsoln[[1]]//Simplify

The function to be minimized is the sum of the product of the gradient squared of each triangle and its area.  The signed area of a triangle (see MathWorld article TriangleArea) is given by:

area[{{x1_,y1_},{x2_,y2_},{x3_,y3_}}]=(1/2)Det[{{x1,y1,1},{x2,y2,1},{x3,y3,1}}]

Therefore the product of the square of the gradient and the area is

grad2Area[{{x1_,y1_,u1_},{x2_,y2_,u2_},{x3_,y3_,u3_}}]=grad2[{{x1,y1,u1},{x2,y2,u2},{x3,y3,u3}}]*Abs[ area[{{x1,y1},{x2,y2},{x3,y3}}] ]//Simplify

Creating the Mesh of Triangles

In this example, the region in question is defined by a list of vertices and the boundary condition is that the value on the boundary is equal to the (((1+y))^2).  The region is shown in the plot p0.

vertices = {{0,0},{(1/2),0},{(1/2),-1},{(5/2),-1},{(5/2),0},{3,0},{3,2},{0,2}};

p0= Graphics[#,AspectRatio -> 1, Axes -> True, AxesOrigin -> {-1,-2}]& @ Line [ vertices /. {a__,b_} -> {a,b,a}]

The first step in creating the mesh is creating the nodes.  Two overlapping arrays of points are used, one from each rectangle.  Union is used to eliminate duplicates.

pts1=Flatten[#,1]& @ Table[{x,y},{x,0,3,(1/4)},{y,0,2,(1/4)}];
pts2=Flatten[#,1]& @ Table[{x,y},{x,(1/2),(5/2),(1/4)},{y,-1,0,(1/4)}];
pts=Union[pts1,pts2];

Show[p0,p1=Graphics[Point /@ pts]]

A function, TriangulatePoints, from the ExtendGraphics package, is used to divide the region into triangles.  The package is available on the Wolfram Library Archive.  As before, triangulation is performed separately for the two rectangles.  This avoids triangles which lie outside the region.

<<ExtendGraphics`Delaunay`

tris1=TriangulatePoints[pts1];
tris2=TriangulatePoints[pts2];
tris=Union[tris1,tris2];

Short[tris]

Graphics[#,AspectRatio -> 1,Axes -> True, AxesOrigin -> {-1,-2}]& @ Join[Line /@ (tris/. {a_,b_,c_}:> {a,b,c,a}),{PointSize[Medium],Point/@ pts}]

Formulating the Problem

It is useful to create list of structures of the form { x, y, u[{ x, y }] } in order to formulate the problem, since grad2area is expecting an argument of this form.

trivals=tris /. {x_?NumberQ,y_?NumberQ} -> {x,y,u[{x,y}]};

The condition that the function u[ { x, y } ] has the value ( y + 1 )^2 on the boundary is written as a rule, which will be used to implement the boundary condition.

bcs=u[{x_,y_}]/; x==0 || x==3 || y==2 || (y==-1 && (1/2)≤x ≤ (5/2)) || (y==0 && (x==0 || x==3)) || (x==(1/2) && -1≤y≤0)||(x==(5/2)&& -1≤y≤0) -> (((y+1))^2);

The problem variables can be extracted from the triangle coordinates using Cases[ ] after the boundary condition has been applied.

Short[vars = Cases[Union @ Flatten[trivals /. bcs],u[_]]]

The value of ((((∂_x) u[x,y]))^2) + ((((∂_y) u[x,y]))^2) is integrated over the region by summing the product of the gradient squared and triangle area.

expr = Plus @@ Simplify[ grad2Area /@ (trivals /. bcs)];

In this situation, it is useful to randomnly perturb the starting values for FindMinimum in order to avoid situations where the gradient of the function to be optimized is zero.

varswithstart=vars/.u[{x_,y_}]->{u[{x,y}],(((y+1))^2)+0.05* RandomReal[]};

The values which minimize the sum are easily found.

soln=FindMinimum[expr,varswithstart,WorkingPrecision -> 16];

Taking the Union of the triangles after applying the solution rule eliminates multiple copies of the points.

res=Union @ Flatten[#,1]& @ (trivals /. bcs /. soln[[2]]);

The results are plotted using two separate ListContourPlots, one for each region, and combining the results.  Two plots are performed in order to avoid getting countour lines outside the region of the two rectangles.

Show[p0,ListContourPlot[res,RegionFunction -> Function[{x,y},#],Contours -> Range[0,20]/2.,ContourShading -> None,ContourLabels -> True]& /@ {(1/2)≤x≤ (5/2) && -1 ≤y ≤ 0,0≤x≤3 &&  0≤y≤2}]

The results compare well to the results found using PDEase2D (shown below), for a calculation with over 12,000 nodes and 6,000 cells.  The largest disagreement is at the contour with value = 1.5, between the two closely spaced corners.