static void
heat_source (int xsize, int ysize, double h[xsize + 2][ysize + 2])
{
  int sizex, sizey;

  // position 4 permanent heat sources
  sizex = xsize / 4;
  sizey = ysize / 4;

  // time-invariant constant values
#pragma omp for parallel
  for (int i = 0; i < sizex; i++)
    for (int j = 0; j < sizey; j++)
      {
    h[1 * xsize / 4 - (sizex / 2) + i][1 * ysize / 4 - (sizey / 2) + j] = HEAT_MAX;
    h[1 * xsize / 4 - (sizex / 2) + i][3 * ysize / 4 - (sizey / 2) + j] = HEAT_MAX;
    h[3 * xsize / 4 - (sizex / 2) + i][1 * ysize / 4 - (sizey / 2) + j] = HEAT_MAX;
    h[3 * xsize / 4 - (sizex / 2) + i][3 * ysize / 4 - (sizey / 2) + j] = HEAT_MAX;
      }
}

/*============================================================================*/
// finite difference method (Laplace)

static void
update (int xsize, int ysize, double h[xsize + 2][ysize + 2], double h_new[xsize + 2][ysize + 2])
{
  // calculate new value based upon neighbor values
#pragma omp parallel
{
#pragma omp for
  for (int i = 1; i <= xsize; i++)
    for (int j = 1; j <= ysize; j++)
      h_new[i][j] = 0.25 * (h[i - 1][j] + h[i + 1][j] + h[i][j - 1] + h[i][j + 1]);

  // update array with new values
#pragma omp for
  for (int i = 1; i <= xsize; i++)
    for (int j = 1; j <= ysize; j++)
      h[i][j] = h_new[i][j];
}
}