ASSIGNMENTS
Assignments for Scientific Computing with DEAL.II

For now, you are free to work on the assignments of your choice, or to propose assignments of interest to you. When turning in a particular assignment, please mention its name, as in "CIRCLE1" or "DEAL.3".

More assignments will be added over time, of course.

• C1: Program in C (Assuming you are not much of a C programmer) Do the LINE1 exercise in C.
• CIRCLE1: Write a program that draws a circle. Display the circle on the computer screen. Print your picture.
• CIRCLE2: Choose 12 points on the circumference of the circle, and let the 13th point be the center. Create 12 triangles by connecting pairs of points on the circumference to the center. Write a program that prints the definition of each of the triangles, by listing the indices of the three points that form it.
• CIRCLE3: Write a program that estimates the area of a circle using your 13 point triangulation from the last assignment, and a formula for the area of a triangle. Compare your answer to the exact answer.
• CIRCLE4: Use the Matlab program called mesh2d to create a mesh for the circle with radius 1. Define the outline of the circle using 12 points. Display your mesh, and print a copy.
• DEAL1: Get deal.ii installed on your computer, or figure out how to use in on the classroom machines.
• DEAL2: Run the deal.ii sample program "step-1"; print the resulting plot.
• DEAL3: If you already got deal.ii set up, use it to define a 2D region shaped like the letter L, define a mesh, and create a plot of the mesh.
• DOC1: Documentation. Get a copy of the "code1.f" program, and a copy of the text book. Look at subroutine PREP, the preprocessor function, and document it, that is, explain the purpose of the subroutine, and how it sets up the problem to be solved. Your documentation should be in the form of comments that are part of the program.
• FORT1: Program in FORTRAN (Assuming you are not much of a FORTRAN programmer) Do the LINE1 exercise in FORTRAN.
• GRAPH1: Graphics project:
  Start with the unit square, and divide it into 4 subsquares, and split each square into two triangles. Throw away one square, so that now you have an "L"-shaped region. Create an "X" file that contains the coordinates of the vertices of the 8 points, and a T file that defines the 6 triangles by listing the indices of triples of points. Write a program that reads in X and T, and plots the region.
• GRAPH2: Graphics project:
  Modify your program so that it can display the index of each point above the point, or the index of each triangle at the center of the triangle. Create a plot with point indices, and a plot with triangle indices, print both and give them to me.
• GRAPH3: Graphics project:
  Modify your program so that each line of the X file includes an X coordinate, a Y coordinate, and a new value we will call U. We might think of U as the temperature at each point. Modify the data file for the L-shaped region by adding U values. You can set each U value in the file to X + Y for now. Read the data into your program. Use some kind of color scheme (such as bluest = coldest, reddest = hottest) to indicate the temperature values. Can you color the points? Can you shade the triangle the average of the three point values? Can you draw contour lines? Can you draw shaded contours?
• GRAPH4: Graphics project:
  Modify the program that you did on tasks 1, 2 and 3 so that it can handle a general region defined by a triangulation. Your program should read data defining points (X,Y) and associated values U, as well as an arrangement of the points into triangles. See if you can display color contours of some kind for the variable U. I can give you some examples of data files.
• LINE1: Begin your 1D finite element program:
  Write "program1", which sets the variables BLEFT and BRIGHT to real numbers, sets N to an integer, defines an array X of size N, fills the array X with evenly spaced values between BLEFT and BRIGHT, and prints them.
• LINE2: Continue your 1D finite element program. Write "program2", which adds new stuff to program1. Set the integer M to N-1. Create a two dimensional array E containing 2 rows and M columns. Each pair of values in E records the indices of the left and right X values for that element. Thus, the first E goes from X(1) to X(2), so E(1,1) = 1 and E(2,1) = 2. The last E goes from X(N-1) to X(N), so E(1,M) = N-1 and E(2,M) = N. Print the values in E. When you have the program set up, try using the same input that program #1 used. Run the program, and save the output.
• LINE3: Continue your 1D finite element program. Write "program3", which adds new stuff to program2. Add a function called BASIS, to evaluate any finite element basis function.

              function value = basis ( t, x, n, j )
              T is a position between BCLEFT and BCRIGHT.
              X is the list of nodes you defined earlier.
              N is the number of nodes.
              J is an index between 1 and N (or 0 and N-1).
            

  Depending on the language you use, the details will differ, but the idea will be roughly that you have to handle the first and last values of J separately; you have to determine which two values of X surround the value of T, and then evaluate the appropriate linear function (up, down, or zero).

              value = 0
              if J = 1, then:
                if T < X(2),            value = ( X(2) - T ) / ( X(2) - X(1) );
              if 1 < J < N then:
                if X(J-1) <  T <= X(J), value = ( T - X(J-1) ) / ( X(J) - X(J-1) ),
                if X(J)   <= T < X(J+1) value = (X(J+1) - T ) / ( X(J+1) - X(J) );
              if J = N then:
                if X(N-1) <  T <= X(N), value = ( T - X(N-1) ) / ( X(N) - X(N-1) ).
            

  Verify that your function is working correctly by having program 3 creating a plot of the values of basis function number J=4, using N = 11 points for the mesh, and 101 equally spaced points T in [BCLEFT,BCRIGHT] for the plot.
• LINE4: Continue your 1D finite element program. Write "program4", which adds new stuff to program 3. Evaluate the interpolant to sin(2*pi*t). We still have BCLEFT=0, BCRIGHT=1, N=11 and our points X. Define a vector C of length N, and set C(J) to be the value of sin(2*pi*X(J)). Our interpolant is

              s(t) = c(1) * b(1)(t) + c(2) * b(2)(t) + ...+ c(n) * b(n)(t).  
            

  where b(1)(t) is the first basis function, evaluated at t, and so on. Compare s(t) and sin(2*pi*t) at 101 equally spaced points in [0,1], with a plot or a table.
• LINE5: Continue your 1D finite element program. Make a copy of program4 and call it program5. We are going to add another tool to the program, which can compute the derivative of any basis function. In program5, make a copy of the basis function code, and call it BASISP. It will look almost the same as basis:

              function value = basisp ( t, x, n, j )
            

  but now we want it to evaluate the derivative of a basis function. Because we are only using piecewise linear functions so far, the derivatives will be very simple to compute:

              if J = 1, then:
                if T < X(2),            value = -1 / ( X(2) - X(1) );
              if 1 < J < N then:
                if X(J-1) <  T <= X(J), value = +1 / ( X(J) - X(J-1) ),
                if X(J)   <= T < X(J+1) value = -1 / ( X(J+1) - X(J) );
              if J = N then:
                if X(N-1) <  T <= X(N), value = +1 / ( X(N) - X(N-1) ).
            

  Verify that your function is working correctly by creating a plot of the derivative of basis function number J=4, using N = 11 points for the mesh, and 101 equally spaced points T in [BCLEFT,BCRIGHT] for the plot.
• LINE6: Numerical integration for 1D. Make a copy of program5 and call it program 6 The linear system that will define our solution involves integrals. We will approximate those integrals using a 2 point quadrature rule:

              integral(a<=x<=b) f(x) dx = (b-a) * ( w1 * f(x1) + w2 * f(x2) ).
            

  Here,
  • w1 = w2 = 1/2;
  • x1 = ((1+s)*a+(1-s)*b)/2;
  • x2 = ((1-s)*a+(1+s)*b)/2;
  • s = 1/sqrt(3).
  Using this rule, and assuming we use N=11 equally spaced nodes on [0,1] as usual, estimate the following integrals, printing your answer each time:
  1. integral ( 0 <= x <= 1 ) cos(2*pi*x) dx;
  3. integral ( ELEMENT #3 ) cos(2*pi*x) dx;
  5. integral ( ELEMENT #3 ) basis(#4) dx;
  7. integral ( ELEMENT #3 ) basisp(#3) dx.
  9. integral ( ELEMENT #3 ) basis(#3) * basis(#4) dx;
  Except for the first integral, we only want to consider the integral over the region that "belongs" to element 3.
• LINE7: The Stiffness Matrix. Make a copy of program6 and call it program 7. Set up a matrix A of size N by N, where N is the number of nodes. We want A(i,j) to be the integral over the whole region of the product of the derivatives of basis functions i and j:

              A(i,j) = integral basisp(i,x) * basisp(j,x) dx
            

  It can be tricky to handle the cases when I=1 and I=N. The cleanest way is to do the following:

              For each ELEMENT
                For I = E(1:2,ELEMENT)
                  For J = E(1:2,ELEMENT)
                    A(i,j) = A(i,j) 
                      + integral ( x in ELEMENT ) basisp(i,x) * basisp(j,x) dx
            

  In other words, do the integral one element at a time. Most of the entries of A(i,j) will be zero. Most of the other entries will be simple, especially if you use N = 11.
• LINE8: The Right Hand Side. Make a new program called program 8, which is a copy of program 7, plus some new stuff. Supposing the right hand side of our differential equation is f(x), we now need to compute a vector called RHS, again using integration. Define a vector RHS of size (N) or, in MATLAB, of size (N,1). We want to compute

              RHS(i) = integral ( 0 <= x <= 1 ) f(x) * basis(i,x) dx
            

  Again, we fill in the values by considering each element:

              For each ELEMENT
                For I = E(1:2,ELEMENT)
                  RHS(i) = RHS(i) 
                    + integral ( x in ELEMENT ) f(x) * basisp(i,x) dx
            

  Use F(X) = -exp(x), use N = 11, compute RHS, and print it.
• LINE9: Boundary conditions. You now have a matrix A and right hand side vector RHS. However, we need to adjust A and RHS to account for boundary conditions at the first node, X(1), and the last node, X(N). All the entries in row 1 of A should be zeroed out, except that A(1,1) should be set to 1. RHS(1) should be set to the boundary value. For this problem, that value is 0. Similarly, set all the entries in row N of the matrix A to zero, except that A(N,N) is 1, and set RHS(N) to 0. Print A and RHS.
• LINE10: Solve the system. Make a new program #10 which is a copy of program #9. We want to solve the linear system

              A*C=RHS
            

  for C, the coefficients of the finite element basis functions. If you are working in MATLAB, this solution step is easy. If you are working in C/C++ or Fortran, you need to locate a linear solver you can use. Print out the 11 entries of your solution vector C.
• LINE11: Tabulate the solution. Copy your previous program and call it program 11. We want to display the solution as a table, and as a graph. Given the coefficients C, the finite element estimate of the solution at any point X is

              UHAT(X) = sum ( 1 <= I <= N ) C(I) * BASIS(I,X)
            

  Use 25 equally spaced points to print a table of UHAT(X). Use 101 equally spaced points to create a plot.
• LINE12: Compare computed and exact solution at nodes. Copy your previous program and call it program 12. Our finite element program was solving the following problem:

              -u'' = - exp(x), u(0) = u(1) = 0.
            

  The exact solution to this problem is

              u(x) = exp(x) + x - x exp(1) - 1
            

  Use 101 equally spaced points to create a "line" plot of your solution, AND a "dot" plot of the exact solution. Expect the line to go right through the dots. The error at the nodes is easy to compute, because at those special places, we have that uhat(x(i)) = c(i), so:

              error = u(x(i)) - uhat(x(i)) = u(x(i)) - c(i)
            

  Print the error at the N=11 nodes.
• LINE13: Compare computed and exact solution as an integral. Copy your previous program and call it program 13. Estimate the following integral:

              L2 = sqrt ( integral ( 0 <= X <= 1 ) ( U(X) - UHAT(X) )^2 dx )
            

  This is a complicated operation. Think about it in small steps, starting with the idea of integrating over a single element.
• LINE14: We want to solve a new problem of the form

              - u'' + 2 * u = 6x^2 - 8x + 4, for 0 <= x <= 1
            

  with boundary conditions

              u(0) = 5, u(1) = 4.
            

  In order to do this, we have to modify our finite element system:

              integral [ u'*b' + 2*u*b ] dx = integral [ (6x^2-8x+4)*b ] dx
            

  We have new boundary conditions and a different right hand side, but the big change is that our matrix A now includes something called the "mass matrix". To set up this new matrix, do the following:

              For each ELEMENT
                For I = E(1:2,ELEMENT)
                  For J = E(1:2,ELEMENT)
                    A(i,j) = A(i,j) 
                      + integral ( x in ELEMENT )     basisp(i,x) * basisp(j,x) dx
                      + integral ( x in ELEMENT ) 2 * basis(i,x)  * basis(j,x)  dx
            

  Using N=11 nodes on [0,1] as usual, set up and solve this new linear system. Compare your computed answer UHAT to the exact answer U(X)=3x^2-4x+5.
• LINE15: Ordinarily, in MATLAB, you solve a linear system A*x=b by writing x=A\b. In other computer languages, you have to write out the instructions to do the solution, or find a prewritten library routine. For the finite element method, one alternative is to use a simple iterative method. The Jacobi iterative method starts with an approximate solution x0 to the linear system, and tries to return an improved estimate x1. The formula is:

              x1(i) = ( b(i) - sum ( j =/= i ) A(i,j) * x0(j) ) / A(i,i)  
            

  Write a function like this: function x1 = jacobi ( A, x0, b ) which does one step of the Jacobi iteration. Let A be the 11x11 matrix which looks like this: [ 2 -1 0 0 0 ...0 0 ] [ -1 2 -1 0 0 ...0 0 ] [ 0 -1 2 -1 0 ...0 0 ] ... [ 0 0 0 0 0 ..-1 2 ] let b be the vector [ 0, 0, 0, ..., 0, 11 ] and let the starting vector x0 be [ 1, 1, 1, ..., 1 ]. The norm of the residual r is the maximum value of |A * x - b|. What is the residual for x0? Compute the improved Jacobi estimate x1. What is the residual for x1? Use x1 as the input to jacobi to get x2, and repeat this process 10 times, so that you get x10. Print the residuals for x0 through x10.
• LINE16: Write a program which starts with the matrix A, right hand side b, and initial guess x0 that we used in the previous example. But now the program must repeatedly call jacobi until the norm of the residual is less than 0.00001. When you satisfy this tolerance, print out the number of iterations you took, and the estimated solution.
• LINE17: Make a copy of PROGRAM14, but now use your Jacobi function to solve the linear system, instead of Matlab's solver. Again, try to get an answer for which the residual is less than 0.00001. How many times did you have to call jacobi?
• MATRIX1: The Wathen matrix is an example of the kinds of matrices that are typically set encountered in finite element calculations. Wathen matrices can be generated of many different sizes, so they are useful when investigating how the accuracy and size limits of a linear solver. Locate a program for generating Wathen matrices. Set up a Wathen matrix A associated with a 10 by 10 spatial grid. How sparse is the matrix? What is its condition? Is it symmetric? Is it positive definite? Try to solve a linear system involving this matrix. (If you do this assignment, then I have some more things for you to try!)
• ODEN1: Copy the file "code1.f" from the class web page. Try to compile it using gfortran. Try to figure out why it will not compile! Ask me for help.
• ODEN2: After we discuss the "code1" program, try harder to get it to compile!
• ODEN3: Try to run the "code1" program using the input file "code1.txt" that is available on the class web page.
• PYTHON1: Program in PYTHON (Assuming you are not much of a Python programmer) Do the LINE1 exercise in PYTHON.
• TRIANGLE1: Begin a 2D finite element program. Create two files that define a triangulation of a region. The "X" file should contain the coordinates of points. The "T" file should contain triples of indices that define the triangles.
• TRIANGLE2: Display your triangulation. Write a program that can read your X and T files, and display the triangulation.
• TRIANGLE3: Read your data into a program. Write a program that reads the X and T files, and which prints out the number of points, the number of triangles, the area of each triangle, and the total area of the region.
• TRIANGLE4: Basis functions. Add a function to your program which can evaluate basis function I at a point (X,Y) in triangle T. If node I is not part of triangle T, the value is zero.

You can return to the DIS 2014 web page.