Page 314 - Numerical Methods for Chemical Engineering

The finite element method 303

or part of an internal boundary between two subdomains. For each of these N e “edge”
line segments in the mesh, the 7×N e array E has a corresponding column. E 1m and E 2m
identify the nodes that “start” and “end” the line segment; their positions are stored in the
corresponding columns of P. We denote a “direction” along the line segment by defining a
contour variable E 3m ≤ s ≤ E 4m such that the segment “starts” at node E 1m with the value
E 3m and ends at node E 2m with the value E 4m . The label E 5m denotes the boundary section
of which this edge segment is a part. Moving in the direction of increasing s, there are
well-defined left-hand and right-hand sides to the line segment. The identity of the region
on the left-hand side (“0” if it is outside of the domain or the subdomain, number 1, 2, . . .
if it is inside) is stored in E 6m and the identity of the region on the right-hand side is stored
in E 7m .
Finally, for each of the N t nonoverlapping triangles in the mesh, there is a corresponding
column in the 4×N t array T. The indices in [1, N p ] of the three nodes that are vertices of
the triangle are stored in T 1m , T 2m , and T 3m . T 4m is the identifier of the subdomain inside
which the triangle lies.
The accuracy of FEM calculations is best for isosceles triangles. To move the node
positions to increase the “quality” of the triangles, type
P = jigglemesh(P, E, T);

Localorglobalmeshrefinementisdoneby refinemesh.Foraglobalrefinement,followed
by quality improvement and plotting, type

[P2, E2, T2] = refinemesh(‘polygon1 geom’, P, E, T);
P2 = jigglemesh(P2, E2, T2); pdemesh(P2, E2, T2);

The algorithms above add points as they deem fit in the interior of the domain using
Delaunay tessellation, an algorithm that can partition a region into non overlapping triangles
for any specified set of nodal positions. For example, from the nodal positions stored in P,
the triangular partition can be constructed from
x = P(1,:)’; y = P(2,:)’;
tri = delaunay(x,y);
triplot(tri,x,y);

Some triangles are formed outside of the domain, but these can be discarded. From the
tessellation, the domain can be partitioned into nonoverlapping Voronoi polyhedra with
nodes at the center of each polyhedron by

voronoi(x,y,tri);
Information about the polyhedra are returned as optional output arguments in voronoi
and voronoin, its extension to three dimensions, four dimensions, etc. (for tessellation, the
extension is delaunayn).
2
2
Finally, various 2-D and 3-D plots of functions such as f (x, y) = x + y can be made
from the nodal function values by pdeplot. pde ex1.m demonstrates the use of this plotting
routine, as well as of the mesh functions described above.