replace node_coordinate and center_coordinate arrays

[alecjohnson]

In the Grid3DCU class, the implementation

  double &getXN(int X, int Y, int Z) { return (node_coordinate[X][Y][Z][0]); }
  double &getYN(int X, int Y, int Z) { return (node_coordinate[X][Y][Z][1]); }
  double &getZN(int X, int Y, int Z) { return (node_coordinate[X][Y][Z][2]); }
  double &getXC(int X, int Y, int Z) { return (center_coordinate[X][Y][Z][0]); }
  double &getYC(int X, int Y, int Z) { return (center_coordinate[X][Y][Z][1]); }
  double &getZC(int X, int Y, int Z) { return (center_coordinate[X][Y][Z][2]); } 

can be replaced with

  const double &getXN(int X, int Y, int Z) { return node_xcoord[X];}
  const double &getYN(int X, int Y, int Z) { return node_ycoord[Y];}
  const double &getZN(int X, int Y, int Z) { return node_zcoord[Z];}
  const double &getXC(int X, int Y, int Z) { return center_xcoord[X];}
  const double &getYC(int X, int Y, int Z) { return center_ycoord[Y];}
  const double &getZC(int X, int Y, int Z) { return center_zcoord[Z];}  

if in the constructor we replace

node_coordinate = newArr4(double, nxn, nyn, nzn, 3);  // 0 -> X, 1 -> Y, 2-> Z
for (int i = 0; i < nxn; i++) {
  for (int j = 0; j < nyn; j++) {
    for (int k = 0; k < nzn; k++) {
      node_coordinate[i][j][k][0] = xStart + (i - 1) * dx;
      node_coordinate[i][j][k][1] = yStart + (j - 1) * dy;
      node_coordinate[i][j][k][2] = zStart + (k - 1) * dz;
    }
  }
}

and

center_coordinate = newArr4(double, nxc, nyc, nzc, 3);
for (int i = 0; i < nxc; i++) {
  for (int j = 0; j < nyc; j++) {
    for (int k = 0; k < nzc; k++) {
      center_coordinate[i][j][k][0] = .5 * (node_coordinate[i][j][k][0] + node_coordinate[i + 1][j][k][0]);
      center_coordinate[i][j][k][1] = .5 * (node_coordinate[i][j][k][1] + node_coordinate[i][j + 1][k][1]);
      center_coordinate[i][j][k][2] = .5 * (node_coordinate[i][j][k][2] + node_coordinate[i][j][k + 1][2]);
    }
  }
}

with

node_xcoord = new double[nxn];
node_ycoord = new double[nyn];
node_zcoord = new double[nzn];
for (int i=0; i<nxn; i++) node_xcoord[i] = xStart + (i - 1) * dx;
for (int j=0; j<nyn; j++) node_ycoord[j] = yStart + (j - 1) * dy;
for (int k=0; k<nzn; k++) node_zcoord[k] = zStart + (k - 1) * dz;
// arrays allocation: cells ---> the first cell has index 1, the last has index ncn-2!
center_xcoord = new double[nxc];
center_ycoord = new double[nyc];
center_zcoord = new double[nzc];
for(int i=0; i<nxc; i++) center_xcoord[i] = .5*(node_xcoord[i]+node_xcoord[i+1]);
for(int j=0; j<nyc; j++) center_ycoord[j] = .5*(node_ycoord[j]+node_ycoord[j+1]);
for(int k=0; k<nzc; k++) center_zcoord[k] = .5*(node_zcoord[k]+node_zcoord[k+1]); 

This reduces memory use and the expense of obtaining coordinate positions.

I surmise that node_coordinate and center_coordinate were created in anticipation of generalizing to a deformed logically cartesian mesh. The rest of the code, however, does not seem to have been consistently written with this in mind. I therefore do not see the value of retaining this implementation, and I would even consider using the alternative accessors

const double &getXN(int X) { return node_xcoord[X];}
const double &getYN(int Y) { return node_ycoord[Y];}
const double &getZN(int Z) { return node_zcoord[Z];}
const double &getXC(int X) { return center_xcoord[X];}
const double &getYC(int Y) { return center_ycoord[Y];}
const double &getZC(int Z) { return center_zcoord[Z];}

based in part on the philosophy that the easiest code to generalize is typically not code that attempts to anticipate generalization but rather code that is clean and simple and readable. For now I will implement these accessors for use in performance-critical parts of the code (pushing particles and summing moments) unless I hear an objection.

[markidis]

Hi Alec,

thanks for this. I agree with you. As you thought, the code was written like this in anticipation of using grid in logical space. thanks, stefano.

[murci3lag0]

Modifications have been included in the amaya-library branch at commit:
5a7a34a