grid.cpp

pic 模拟程序！面向对象

	for (j = 0; j < J; j++)
	{
		xlo = X[j][0].e1(); xhi = X[j+1][0].e1();
		if((x1 >= xlo) && (x1 <= xhi)) break;
	}
	dx=xhi-xlo;
	if(j==J) xlo = X[j][0].e1(); // catch things near the last edge.
	x1grid = j + (x1 - xlo)/dx;

	// round it off to the nearest grid
	if(x1grid >= 0) x1grid = x1grid + 0.5;
	else x1grid = x1grid - 0.5;

	if(x1grid > J+0.75  || x1grid < -0.75) { x1grid = -1; }  // outside this mesh.

	Scalar x2 = xMKS.e2(), x2grid;
	int k;
	for (k = 0; k < K; k++)
	{
		xlo = X[0][k].e2(); xhi = X[0][k+1].e2();
		if((xlo <= x2) && (x2 <= xhi)) break;
	}
	dx=xhi-xlo;
	if(k==K) xlo = X[0][k].e2(); // catch things near the last edge.
	x2grid = k + (x2 - xlo)/dx;


	// round it off to the nearest grid
	if(x2grid >= 0) x2grid = x2grid + 0.5;
	else x2grid = x2grid - 0.5;

	if(x2grid > K +0.75  || x2grid < -0.75) { x2grid = -1; }  // outside this mesh.

	if(x1grid == -1 || x2grid == -1) return Vector2(-1,-1);
	// truncate
	x1grid = (int) x1grid;
	x2grid = (int) x2grid;
	return Vector2(x1grid, x2grid);						//	replace with Westermann method.
}

//-----------------------------------------------
//  Initialize the cell volumes array.  This array is
//  used by the charge collection routine.

void Grid::initHalfCellVolumes(void)
{
	int	j, k;

	switch(geometryFlag)
	{
		case ZXGEOM:
			for (j=0; j<=J; j++)
			{
				Scalar	zp, zm, xp;
				Scalar	xm = X[j][0].e2();
				for (k=0; k<=K; k++)
				{
					zp = (j==J) ? X[j][k].e1() : X[j+1][k].e1();
					zm = (j==0) ? X[j][k].e1() : X[j-1][k].e1();
					xp = (k==K) ? X[j][k].e2() : X[j][k+1].e2();
					xm = (k==0) ? X[j][k].e2() : X[j][k-1].e2();
					halfCellVolumes[j][k] = 0.25*(zp - zm)*(xp - xm);
				}
			}
		break;

		case ZRGEOM:
		default:
			for (j=0; j<=J; j++)
			{
				Scalar	zp, zm, rp2;
				Scalar	rm2 = sqr(X[j][0].e2());
				for (k=0; k<=K; k++)
				{
					zp = (j==J) ? X[j][k].e1() : X[j+1][k].e1();
					zm = (j==0) ? X[j][k].e1() : X[j-1][k].e1();
					rp2 = (k==K) ? sqr(X[j][k].e2())
						: sqr(0.5*(X[j][k+1].e2() + X[j][k].e2()));
					halfCellVolumes[j][k] = M_PI*0.5*(zp - zm)*(rp2 - rm2);
					if (axis()) // apply radial half cell vol corrections
					  if (k == 0) halfCellVolumes[j][k] /= radial0;
					  else if (k==K) halfCellVolumes[j][k] /= radialK;
					rm2 = rp2;
				}
			}
		break;
		}
	if (PeriodicFlagX1)
		for (k=0;k<=K;k++){
			halfCellVolumes[0][k] += halfCellVolumes[J][k];
			halfCellVolumes[J][k] = halfCellVolumes[0][k];
		}
	if (PeriodicFlagX2)
		for (j=0; j<=J; j++){
			halfCellVolumes[j][0] += halfCellVolumes[j][K];
			halfCellVolumes[j][K] = halfCellVolumes[j][0];
		}
}

void Grid::init_cellVolume(void)
{
	int j,k;
	switch(geometryFlag)
		{
		case ZXGEOM: 
			for (j=0; j<J; j++)
				for (k=0; k<K; k++)
					cellVolume_array[j][k] = (X[j][k+1].e2() - X[j][k].e2())*
						(X[j+1][k].e1() - X[j][k].e1());
			break;
		case ZRGEOM:
		default: 
			for (j=0; j<J; j++)
				for (k=0; k<K; k++)
					cellVolume_array[j][k] = M_PI*(X[j][k+1].e2()*X[j][k+1].e2()
						- X[j][k].e2()*X[j][k].e2())*(X[j+1][k].e1()	- X[j][k].e1());
			break;
		}
}

void Grid::init_dl(void)
{
	int j,k;
	switch(geometryFlag)
		{
		case ZXGEOM: 
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dl_array[j][k] = Vector3(((j==J) ? 0 : X[j+1][k].e1() - X[j][k].e1()),
																	 ((k==K) ? 0 : X[j][k+1].e2() - X[j][k].e2()),
																	 1.0);
					
			break;
		case ZRGEOM:
		default: 
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dl_array[j][k] = Vector3(((j==J) ? 0 : X[j+1][k].e1() - X[j][k].e1()),
														 ((k==K) ? 0 : X[j][k+1].e2() - X[j][k].e2()),
														 2*M_PI*X[j][k].e2());
			break;
		}
	if (PeriodicFlagX1)
		for (k=0;k<=K;k++)
			dl_array[J][k] = dl_array[0][k];
		
	if (PeriodicFlagX2)
		for (j=0; j<=J; j++)
			dl_array[j][K] = dl_array[j][0];
}

void Grid::init_dS(void)
{
	int j,k;
	switch(geometryFlag)
		{
		case ZXGEOM: 
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dS_array[j][k] = Vector3(((k==K) ? 0 : X[j][k+1].e2() - X[j][k].e2()),
														 ((j==J) ? 0 : X[j+1][k].e1()-X[j][k].e1()),
														 (((j==J)||(k==K)) ? 0 : (X[j+1][k].e1() - X[j][k].e1())*(X[j][k+1].e2() - X[j][k].e2())));
			break;
		case ZRGEOM:
		default: 
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dS_array[j][k] = Vector3(((k==K) ? 0 : M_PI*(X[j][k+1].e2()*X[j][k+1].e2() - X[j][k].e2()*X[j][k].e2())),
														 ((j==J) ? 0 : 2*M_PI*X[j][k].e2()*(X[j+1][k].e1() - X[j][k].e1())),
														 (((j==J)||(k==K)) ? 0 : (X[j+1][k].e1() - X[j][k].e1())*(X[j][k+1].e2() - X[j][k].e2())));
			break;
		}
	if (PeriodicFlagX1)
		for (k=0;k<=K;k++)
			dS_array[J][k] = dS_array[0][k];
		
	if (PeriodicFlagX2)
		for (j=0; j<=J; j++)
			dS_array[j][K] = dS_array[j][0];
}

void Grid::init_dlPrime(void)
{
	int j,k;
//	Scalar dlPrime1, dlPrime2;
	switch(geometryFlag)
		{
		case ZXGEOM: 
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dlPrime_array[j][k] = Vector3(0.5*(((j==J) ? X[j][k].e1() : X[j+1][k].e1()) 
																						 - ((j==0) ? X[j][k].e1() : X[j-1][k].e1())),
																				0.5*(((k==K) ? X[j][k].e2() : X[j][k+1].e2()) 
																						 - ((k==0) ? X[j][k].e2() : X[j][k-1].e2())),
																				1.0);

			break;
		case ZRGEOM:
		default: 
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dlPrime_array[j][k] = Vector3(0.5*(((j==J) ? X[j][k].e1() : X[j+1][k].e1()) 
																						 - ((j==0) ? X[j][k].e1() : X[j-1][k].e1())),
																				0.5*(((k==K) ? X[j][k].e2() : X[j][k+1].e2()) 
																						 - ((k==0) ? X[j][k].e2() : X[j][k-1].e2())),
																				(k==K) ? 2*M_PI*X[j][k].e2() : M_PI*(X[j][k+1].e2() + X[j][k].e2()));
			break;
		}
	if (PeriodicFlagX1)
		for (k=0;k<=K;k++){
			dlPrime_array[0][k].set_e1(dlPrime_array[0][k].e1()+dlPrime_array[J][k].e1());
			dlPrime_array[J][k].set_e1(dlPrime_array[0][k].e1());
		}
	if (PeriodicFlagX2)
		for (j=0; j<=J; j++){
			dlPrime_array[j][0].set_e2(dlPrime_array[j][0].e2()+dlPrime_array[j][K].e2());
			dlPrime_array[j][K].set_e2(dlPrime_array[j][0].e2());
		}

}

void Grid::init_dSPrime(void)
{
	int j,k;
	// all geometry switches are handled by dlPrime
	switch(geometryFlag)
		{
		case ZXGEOM:
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dSPrime_array[j][k] = Vector3(dlPrime_array[j][k].e2()*dlPrime_array[j][k].e3(),
																				dlPrime_array[j][k].e1()*dlPrime_array[j][k].e3(),
																				dlPrime_array[j][k].e1()*dlPrime_array[j][k].e2());
			break;
		case ZRGEOM:
			for (j=0; j<=J; j++)
				for (k=0; k<=K; k++)
					dSPrime_array[j][k] = Vector3((k==K) ? M_PI*(sqr(X[j][k].e2())-sqr(.5*(X[j][k].e2() + X[j][k-1].e2()))) : 
																				(k==0) ? M_PI*(sqr(.5*(X[j][k+1].e2() + X[j][k].e2()))-sqr(.5*(X[j][k].e2()))) : 
																				M_PI*(sqr(.5*(X[j][k+1].e2() + X[j][k].e2()))-sqr(.5*(X[j][k].e2() + X[j][k-1].e2()))),
																				dlPrime_array[j][k].e1()*dlPrime_array[j][k].e3(),
																				dlPrime_array[j][k].e1()*dlPrime_array[j][k].e2());
			break;
		}
// Periodic handled by dlPrime_array
}


#ifdef MPI_VERSION

extern int MPI_RANK;
int Grid::getGlobalIndex(int j,int k, GRIDDIR dir) {
  // first compute the index of the j,k node
  int ROW = MPI_RANK/M, COLUMN=MPI_RANK%M;
  int rejstart = MPIgstarts[COLUMN][ROW];
  int result = rejstart + j + k*(J+1);
  int test;

  switch(dir) {
  case HERE: 
	 {
	 	return result;
	 }
  case UP:  // lesser k, i.e., above.
	 {
		test = result - (J+1);
		if(test >= rejstart) return test;
		// OK, the point above us is outside this region.
		// if we're in the top row, return an error condition
		if(ROW==0) return -1;
		return (int)( MPIgstarts[COLUMN][ROW -1] + 
		  (MPInodeDim[COLUMN][ROW-1].e1()+1)*(MPInodeDim[COLUMN][ROW-1].e2()-1) + j);
	 }
  case RIGHT:  //greater j, i.e., right.
	  {
		  if(j<J) return result+1;  // it's inside the system
		  if(M-COLUMN == 1) return -1;  // it's outside the simulation!!
		  return (int)(MPIgstarts[COLUMN+1][ROW] + k * (MPInodeDim[COLUMN+1][ROW].e1()+1) +1);
	  }
  case DOWN: // greater k, i.e., down.
	  {
		  if(k<K) return result + (J+1);
	  if(N-ROW == 1) return -1;  //it's outside the simulation!!
		  return (int)(MPIgstarts[COLUMN][ROW+1] + (MPInodeDim[COLUMN][ROW+1].e1()+1) +j);
	  }
  case LEFT: //lesser j, i.e., left.
	  {
		  if(j>0) return result -1;
		  if(COLUMN==0) return -1;
		  return (int)(MPIgstarts[COLUMN-1][ROW] + (k+1)*(MPInodeDim[COLUMN-1][ROW].e1()+1)-2);
	  }
  case UPWRAP:  // wrapping around.
	  {
		  return (int)(MPIgstarts[COLUMN][N-1] + 
			  (MPInodeDim[COLUMN][N-1].e1()+1)*(MPInodeDim[COLUMN][N-1].e2()-1) +j);
	  }
  case RIGHTWRAP: // wrapping around
	  {
		 return (int)(MPIgstarts[0][ROW] + k * (MPInodeDim[0][ROW].e1()+1) +1);
	 } 
  case DOWNWRAP: // wrapping around
	  {
		return (int)(MPIgstarts[COLUMN][0] + (MPInodeDim[COLUMN][0].e1()+1) +j);  
	}
  case LEFTWRAP: // wrapping around.
	  {
		  return (int)(MPIgstarts[M-1][ROW] + (k+1)*(MPInodeDim[M-1][ROW].e1()+1)-2);
	  }
  default:
  return result;
  }
}

#endif //MPI_VERSION