The data structure is described in the M-file Mesh1.m (this M-file contains only comments, which are displayed when you type "help Mesh1"):
addpath('fempack')
help Mesh1
The Fem1 package implements piecewise linear finite elements for two
dimensional problems and is intended to accompany "Partial
Differential Equations: Analytical and Numerical Methods" (second
edition) by Mark S. Gockenbach (SIAM 2010). It uses the data
structure described in Section 13.1 of the text (that is, the four
arrays NodeList, NodePtrs, FNodePtrs, and ElList, augmented by
three more arrays (CNodePtrs, ElEdgeList, and FBndyList) to allow
the code to handle inhomogeneous boundary conditions (as hinted at
in the text). These seven arrays are collected into a structure.
NodeList: Mx2 matrix, where M is the number of nodes in the mesh
(including boundary nodes). Each row corresponds to a
node and contains the coordinates of the node.
NodePtrs: Mx1 matrix. Each entry corresponds to a node:
if node i is free, NodePtrs(i) is the index of node i
in FNodePtrs; else NodePtrs(i) is the negative of
index of node i in CNodePtrs
FNodePtrs: Nx1 vector, where N is the number of free nodes.
NodePtrs(i) is the index into NodeList of the ith
free node. So, for example, NodeList(NodePtrs(i),1:2)
are the coordinates of the ith free node.
CNodePtrs: Kx1 vector, where K=M-N is the number of constrained
nodes. CNodePtrs(i) is the index into NodeList of the
ith constrained node.
ElList: Lx3 matrix, where L is the number of triangular elements
in the mesh. Each row corresponds to one element and
contains pointers to the nodes of the triangle in
NodeList.
ElEdgeList: Lx3 matrix. Each row contains flags indicating
whether each edge of the corresponding element is on
the boundary (flag is -1 if the edge is a constrained
boundary edge, otherwise it equals the index of the
edge in FBndyList) or not (flag is 0). The edge
of the triangle are, in order, those joining vertices
1 and 2, 2 and 3, and 3 and 1.
FBndyList: Bx2 matrix, where B is the number of boundary edges
not constrained by Dirichlet conditions. Each
row corresponds to one edge and contains pointers
into NodeList, yielding the two vertices of the edge.
The command RectangleMeshD1 creates a regular triangulation of a rectangle
,
assuming Dirichlet conditions on the boundary.
help RectangleMeshD1
T=RectangleMeshD1(nx,ny,lx,ly)
This function creates a regular, nx by ny finite element mesh for
a Dirichlet problem on the rectangle [0,lx]x[0,ly].
The last three arguments can be omitted; their default values are
ny=nx, lx=1, ly=1. Thus, the command "T=RectangleMeshD1(m)"
creates a regular mesh with 2m^2 triangles on the unit square.
For a description of the data structure describing T, see
"help Mesh1".
The command RectangleMeshN1 creates a regular triangulation of the same rectangle, assuming Neumann conditions on the boundary, while RectangleMeshTopD1 imposes Dirichlet conditions on the top edge and Neumann conditions on the other three edges.
Thus I only provide the means to deal with a single domain shape (a rectangle), and only under three combinations of boundary conditions. To use my code to solve BVPs on another domain, you will have to write code to generate the mesh yourself. Note, however, that the rest of the code is general: it handles any polygonal domain, with any combination of homogeneous Dirichlet and Neumann conditions.
Here I create a mesh:
clear
T=RectangleMeshD1(4)
T =
NodeList: [25x2 double]
NodePtrs: [25x1 double]
FNodePtrs: [9x1 double]
CNodePtrs: [16x1 double]
ElList: [32x3 double]
ElEdgeList: [32x3 double]
FBndyList: []
The mesh I just created is shown in Figure 10.1 in the text. I have provided a means of viewing a mesh:
help ShowMesh1
ShowMesh1(T,flag)
This function displays a triangular mesh T. Unless the flag
is nonzero, free nodes are indicated by a 'o', and
constrained nodes by an 'x'.
For a description of the data structure T, see "help Mesh1".
The optional argument flag has the following effect:
flag=1: the triangles are labeled by their indices
flag=2: the nodes by their indices
flag=3: both the nodes and triangles are labeled
flag=4: the free nodes are labeled by their indices.
flag=5: the free and constrained nodes are labeled by
their indices.
ShowMesh1(T)
(Notice that it is also possible to label the triangles and/or nodes by their indices.)
Share with your friends: |