Assignment 3

Consider a rectangular piece of metal alloy, two times as wide as high, consisting of three different metals, each with different thermal characteristics. For each region of the alloy, there is a given amount (expressed in terms of a percentage) of each of the three base metals, that varies up to 25 percent due to random noise. The top left corner (at the mesh element at index [0,0]) is heated at $S$ degrees Celsius and the bottom right corner (index [width - 1,height - 1]) is heated at $T$ degrees Celsius. The temperature at these points may also randomly vary over time.

Your program calculates the final temperature for each region on the piece of alloy. The new temperature for a given region of the alloy is calculated using the formula

\begin{displaymath}temp = \sum_{m = 1}^{3} C_m * \left(\sum_{n \in N} temp_n *
p^{m}_{n}\right) / \vert N\vert\end{displaymath}

where $m$ represents each of the three base metals, $C_m$ is the thermal constant for metal $m$, $N$ is the set representing the neighbouring regions, $temp_n$ is the temperature of the neighbouring region, $p^{m}_{n}$ is the percentage of metal $m$ in neighbour $n$, and $\vert N\vert$ is the number of neighbouring regions. This computation is repeated until the temperatures converge to a final value or a reasonable maximum number of iterations is reached.

The values for $S$, $T$, $C_1$, $C_2$, $C_3$, the dimensions of the mesh, and the threshold should be parameters to the program. Note however, that combinations of these parameters do not do not converge well. Try values of (0.75, 1.0, 1.25) C1, C2, C3 for your test/demo.

Assume that the edges are maximally insulated from their surroundings.

Your program should graphically display the results by drawing a grid of points with intensities or colors indicating temperature. (Alternatively, you can just output them to a file and use something like gnuplot to display them.)

Acknowledgment: This assignment was adapted from a study by Steve MacDonald.


Here's is some code for the plain Jacobi example in CPJ:
abstract class JTree extends RecursiveAction {
  volatile double maxDiff; // for convergence check

class Interior extends JTree {
  private final JTree[] quads;

  Interior(JTree q1, JTree q2, JTree q3, JTree q4) {
    quads = new JTree[] { q1, q2, q3, q4 };

  public void compute() { 
    double md = 0.0;
    for (int i = 0; i < quads.length; ++i) {
      md = Math.max(md,quads[i].maxDiff);
    maxDiff = md;

class Leaf extends JTree {
  private final double[][] A; private final double[][] B; 
  private final int loRow;    private final int hiRow;
  private final int loCol;    private final int hiCol; 
  private int steps = 0; 

  Leaf(double[][] A, double[][] B, 
       int loRow, int hiRow, int loCol, int hiCol) {
    this.A = A;   this.B = B;
    this.loRow = loRow; this.hiRow = hiRow;
    this.loCol = loCol; this.hiCol = hiCol;

  public void compute() {
    boolean AtoB = (steps++ % 2) == 0;
    double[][] a = (AtoB)? A : B;
    double[][] b = (AtoB)? B : A;
    double md = 0.0;
    for (int i = loRow; i <= hiRow; ++i) {
      for (int j = loCol; j <= hiCol; ++j) {
        b[i][j] = 0.25 * (a[i-1][j] + a[i][j-1] +
                          a[i+1][j] + a[i][j+1]);
        md = Math.max(md, Math.abs(b[i][j] - a[i][j]));
    maxDiff = md;

class Jacobi extends RecursiveAction {
  static final double EPSILON = 0.001; // convergence criterion
  final JTree root; 
  final int maxSteps;
  Jacobi(double[][] A, double[][] B, 
         int firstRow, int lastRow, int firstCol, int lastCol,
         int maxSteps, int leafCells) {
    this.maxSteps = maxSteps;
    root = build(A, B, firstRow, lastRow, firstCol, lastCol,

  public void compute() {
    for (int i = 0; i < maxSteps; ++i) {
      if (root.maxDiff < EPSILON) { 
      else root.reinitialize();

  static JTree build(double[][] a, double[][] b,
                     int lr, int hr, int lc, int hc, int size) {
    if ((hr - lr + 1) *  (hc - lc + 1) <= size) 
      return new Leaf(a, b, lr, hr, lc, hc);
    int mr = (lr + hr) / 2; // midpoints
    int mc = (lc + hc) / 2;
    return new Interior(build(a, b, lr,   mr, lc,   mc, size),
                        build(a, b, lr,   mr, mc+1, hc, size),
                        build(a, b, mr+1, hr, lc,   mc, size),
                        build(a, b, mr+1, hr, mc+1, hc, size));