22M:176 Finite Elements
Programming Assignment (Assignment #3)
Due Wednesday, December 9th, 1998

In this assignment we will build a simple FEM solver for two-dimensional PDE's of the form

$$\begin{array}{cl} \nabla\cdot a({\bf x})\nabla u + c({\bf ... ...f x}) = g({\bf x}),& \qquad {\bf x}\in\partial\Omega\end{array}$$

We will use MATLAB to solve the system of equations, and the triangulation will be provided.

We will use piecewise linear basis functions under the given triangulation. The basis function $\phi_i$ is the piecewise linear function which is one at vertex ${\bf v}_i$ of the triangulation, and zero at all other vertices of the triangulation.

So, our purpose is to define the system of linear equations to solve. We will need to compute the integrals

$$\begin{array}{c}\displaystyle \int_\Omega[a\nabla\phi_i\c... ...phi_j] dV,\\ \displaystyle \int_\Omega f \phi_i dV.\end{array}$$

To approximate these integrals we will use the following quadrature formula for a triangle K with vertices ${\bf v}_i$, i=1,2,3:

$$\int_K h dV \approx \vert K\vert\left[ \frac1{20}\sum_i h({... ...um_{i,j:i<j}h({\bf v}_{ij}) +\frac9{20}h({\bf v}_{123})\right]$$

where ${\bf v}_{ij}=({\bf v}_{i}+{\bf v}_j)/2$,${\bf v}_{123}=({\bf v}_1+{\bf v}_2+{\bf v}_3)/3$, and |K| is the area of triangle K.

Data structures

Perhaps the most complex part of the coding is organizing the data in an efficient way. MATLAB has sparse matrices, which you should use for storing the linear equations.

The triangulation will be represented by a collection of arrays as follows: Suppose that there are N nodes (including both internal and boundary nodes) and M triangles.

• vertex: This will be an $N\times2$ array of real numbers where vertex(i,1) is the x-component of ${\bf v}_i$, and vertex(i,2) is its y-component.
• triangle: This will be an $M\times3$ array of integers where triangle(j,1), triangle(j,2), and triangle(j,3) are the indexes of the vertices of the triangle K_j.
• interior: This is an N component array of integers; interior(i) is 1 if ${\bf v}_i$ is an interior node and zero if it is a boundary node.
• triang1 and next_triang: These two arrays give an efficient way of determining which triangles touch a given vertex. triang1 is an N component array where triang1(i) is the index (in next_triang) of the ``first'' triangle touching vertex ${\bf v}_i$. next_triang is a two-column array which translates the entries in triang1 to indexes into triangle, as well as giving the index (into next_triang) of the following triangles touching ${\bf v}_i$. Here is some MATLAB code to process all triangles touching vertex ${\bf v}_i$:

  % Find all triangles touching vertex # i
  tri_index = triang1(i);
  while tri_index > 0
    k = next_triang(tri_index,1);
    % Process triangle # k
      ..... triangle(k,1) ....
    tri_index = next_triang(tri_index,2);
  end

These arrays will be computed by a routine mk_nxt_tri.

For computing integrals involving $\phi_i$, we need to find all of the triangles K that touch ${\bf v}_i$.This can be done using the above code. For computing

$$\int_\Omega[a\nabla\phi_i\cdot\nabla\phi_j+c\phi_i\phi_j] dV$$

we should only integrate over only the triangles that touch both${\bf v}_i$ and ${\bf v}_j$.Since each vertex is touched by only six or so triangles, we can test each triangle as we find it using the above code:

  % Find all triangles touching vertex # i and vertex # j
  tri_index = triang1(i);
  while tri_index > 0
    k = next_triang(tri_index,1);
    if triangle(k,1) == j | triangle(k,2) == j | triangle(k,3) == j
      % Integrate over triangle # k
        ..... triangle(k,1) ....
    end
    tri_index = next_triang(tri_index,2);
  end

How do we find |K|?

The area of a triangle K with vertices ${\bf v}_1$, ${\bf v}_2$ and${\bf v}_3$ can be computed as$\frac12\vert({\bf v}_1-{\bf v}_3)\times({\bf v}_2-{\bf v}_3)\vert$ where ${\bf a}\times{\bf b}$ is the vector cross product. Alternatively, if ${\bf a}=(a_x,a_y)$ in co-ordinate form,

$$\vert K\vert = \frac12\vert(v_{1x}-v_{3x})(v_{2y}-v_{3y})-((v_{1y}-v_{3y})(v_{2x}-v_{3x})\vert.$$

since ${\bf a}\times{\bf b}= (a_xb_y-a_yb_x){\bf k}$ where ${\bf k}$ is the unit vector in the ``z'' direction.

Boundary conditions

When you formulate the Galerkin equations to solve, write

$$U({\bf x}) = \sum_{i\in{\cal I}}u_i\phi_i({\bf x}) + \sum_{i\in{\cal B}}u_i\phi_i({\bf x})$$

where $\cal I$ is the set of interior nodes and $\cal B$ is the set of boundary nodes. If $i\in{\cal B}$, make the i'th equation just $u_i = g({\bf v}_i)$.For $i\in{\cal I}$, the i'th equation should be the usual equation:

$$\sum_{j=1}^N c_{ij}u_j = \int_\Omega f \phi_i dV.$$

Test problem

We will have two test problems. The first will be to solve the Poisson equation $-\nabla^2u=1$ on the unit square $\Omega=(0,1)\times(0,1)$ with homogeneous Dirichlet boundary conditions (u=0 on $\partial\Omega$). For comparison purposes, the exact solution is

$$u(x,y) = \sum_{i,j=0}^\infty a_{ij}\sin(i\pi x)\sin(j\pi y)$$

for suitable coefficients a_ij. Find the correct coefficients for this problem.

The triangulation to use is the regular triangulation described in Q.3 or the 2nd assignment. Use 11 grid points on each axis.

The second test problem will be defined by the data in the ASCII files vertex, triangles, triang1, and next_tri which contain the data for the corresponding arrays. These are available via the WWW page . The mesh described by the above files looks like this:

ap_mesh.pdf

11/18/1998