C++ Grid3D类代码示例

void MeshFileManager::writeShape(const Grid3D &grid, const QString &path)
{
    QFile file(getRootPath() + path);
    if (!file.open(QIODevice::WriteOnly)) {
        return;
    }

    QTextStream stream(&file);
    stream  << "MeshVersionFormatted 1" << endl
            << "Dimension" << endl
            << "3" << endl;

    stream  << "Vertices" << endl
            <<  grid.getNx() * grid.getNy() * grid.getNz() << endl;

    for (int z = 0; z < grid.getNz(); ++z) {
        for (int y = 0; y < grid.getNy(); ++y) {
            for (int x = 0; x < grid.getNx(); ++x) {
                Point3D currentPoint(grid.getOrigin().getX() + x,
                                     grid.getOrigin().getY() + y,
                                     grid.getOrigin().getZ() + z);

                if (grid.getProperty(currentPoint.getX(), currentPoint.getY(), currentPoint.getZ()) < 0) {
                    stream << currentPoint.format()
                           << endl;
                }
            }
        }
    }

    stream << "End" << endl;
}

void MeshFileManager::writeGeometry(const Grid3D &grid, const QString &path)
{
    QFile file(getRootPath() + path);
    if (!file.open(QIODevice::WriteOnly)) {
        return;
    }

    QTextStream stream(&file);
    stream  << "MeshVersionFormatted 1" << endl
            << "Dimension" << endl
            << "3" << endl;

    int nx = grid.getNx();
    int ny = grid.getNy();
    int nz = grid.getNz();
    Point3D origin = grid.getOrigin();

    stream  << "Vertices" << endl
            <<  nx * ny * nz << endl;
    QVector<Hexaedra> hexs;

    int xD = 1;
    int yD = nx;
    int zD = nx * ny;
    for (int z = 0; z < nz; ++z) {
        for (int y = 0; y < ny; ++y) {
            for (int x = 0; x < nx; ++x) {
                stream << Point3D(origin.getX() + x, origin.getY() + y, origin.getZ() + z).format()
                       << endl;

                if (!(x == nx - 1 || y == ny - 1 || z == nz - 1)) {
                    hexs.push_back(Hexaedra(
                        xD * (x + 0) + yD * (y + 0) + zD * (z + 0),
                        xD * (x + 1) + yD * (y + 0) + zD * (z + 0),
                        xD * (x + 1) + yD * (y + 1) + zD * (z + 0),
                        xD * (x + 0) + yD * (y + 1) + zD * (z + 0),
                        xD * (x + 0) + yD * (y + 0) + zD * (z + 1),
                        xD * (x + 1) + yD * (y + 0) + zD * (z + 1),
                        xD * (x + 1) + yD * (y + 1) + zD * (z + 1),
                        xD * (x + 0) + yD * (y + 1) + zD * (z + 1)
                    ));
                }
            }
        }
    }

    stream  << "Hexahedra" << endl
            << hexs.size() << endl;

    for (int i = 0; i < hexs.size(); ++i) {
        stream << hexs.at(i).format() << endl;
    }

    stream << "End" << endl;
}

void Pic::kick(Grid3D& Ex, Grid3D& Ey, double ds){
  double beta0 = SP->beta0, gamma0 = SP->gamma0, A = SP->A, Z = SP->Z;
  double temp = qe*Z/(mp*A*pow(gamma0, 3)*pow(beta0*clight, 2))*ds;
  double ex, ey;

  for(long j=0; j<pics.size(); ++j){
    ex = temp*Ex.Grid2PIC(pics[j].x, pics[j].y, pics[j].z);
    ey = temp*Ey.Grid2PIC(pics[j].x, pics[j].y, pics[j].z);

  //momenta:

    pics[j].xs += ex;
    pics[j].ys += ey;
  }
}

void poisson_xyz(Grid3D& Ex, Grid3D& Ey, Grid3D& rho, Greenfb& gf)
{
 int NX=rho.get_NX(),NY=rho.get_NX(), NZ=rho.get_NZ();
 double dx=rho.get_dx(), dy=rho.get_dy();
 for(int i=0; i<NZ; i++)
 {
  poisson_xy(rho[i],gf);  
  for(int j=1;j<NX-1;j++)
   for(int l=1;l<NY-1;l++)
    {
     Ex[i](j,l)=-(rho[i](j+1,l)-rho[i](j-1,l))/(2.0*dx);
     Ey[i](j,l)=-(rho[i](j,l+1)-rho[i](j,l-1))/(2.0*dy);
    }
 }
}

/*
typedef GridSerializerVTK< ls::Grid3D<double>, StructuredPointsTraits > BasicSerializerVTK_SPoint;
// The list of types we want to test.
typedef testing::Types<BasicSerializerVTK, BasicSerializerVTK_SPoint> Implementations;

TYPED_TEST_CASE(VTKTest, BasicSerializerVTK);
*/
TEST(VTKTest, writeAndRead3)
{
  Close_relative close_at_tol(10e-8); // for whatever reason vtk writes with error

  size_t n = 12, m = 11, w = 14;
  Grid3D<double> grid(n, m, w);

  double top[] = {4.0, 5.0, 9.0};
  double low[] = {-3.0, -4.0, -5.0};
  Box3D box(low, top);
  double h[] = {box.getSizeX() / (n - 1.0), box.getSizeY() / (m - 1.0), box.getSizeZ() / (w - 1.0)};

  for (size_t i = 0; i < n; ++i) {
    for (size_t j = 0; j < m; ++j) {
      for (size_t k = 0; k < w; ++k) {
        MathVector3D p(i * h[0], j * h[1], k * h[2]);
        p += box.getLow();
        assert(box.inside(p));
        grid(i, j, k) = testFunction2(p);
      }
    }
  }

  grid.setBoundingBox(box);
  BasicSerializerVTK writer(grid, "test-aux/writeAndRead3");
  writer.run();

  Grid3D<double> readGrid; // at this point I don't know the size
  BasicDeserializerVTK reader(readGrid, "test-aux/writeAndRead3");
  reader.run();

  EXPECT_EQ(grid.getBoundingBox(), readGrid.getBoundingBox());

  for (size_t i = 0; i < n; ++i) {
    for (size_t j = 0; j < m; ++j) {
      for (size_t k = 0; k < w; ++k) {
        if (grid(i, j, k) != readGrid(i, j, k)) {
          EXPECT_TRUE(close_at_tol(grid(i, j, k), readGrid(i, j, k)));
        }
      }
    }
  }
}

TEST(VTKTest, writeAndRead1)
{
  Close_relative close_at_tol(10e-8); // for whatever reason vtk writes with error

  size_t n = 8, m = 8, w = 8;
  Box3D domainWrite(1.0, 1.0, 1.0);
  Grid3D<double> grid(n, m, w);
  double h = 1.0 / (n - 1);
  for (size_t i = 0; i < n; ++i) {
    for (size_t j = 0; j < m; ++j) {
      for (size_t k = 0; k < w; ++k) {
        double point[] = { i * h - 0.5, j * h - 0.5, k * h - 0.5};
        grid(i, j, k) = testFunction1(point);
      }
    }
  }

  grid.setBoundingBox(domainWrite);
  BasicSerializerVTK writer(grid, "test-aux/writeAndRead1");
  EXPECT_TRUE(writer.run());

  Grid3D<double> readGrid; // at this point I don't know the size
  BasicDeserializerVTK reader(readGrid, "test-aux/writeAndRead1");
  EXPECT_TRUE(reader.run());

  EXPECT_EQ(grid.getBoundingBox(), readGrid.getBoundingBox());

  for (size_t i = 0; i < n; ++i) {
    for (size_t j = 0; j < m; ++j) {
      for (size_t k = 0; k < w; ++k) {
        if (grid(i, j, k) != readGrid(i, j, k)) {
          EXPECT_TRUE(close_at_tol(grid(i, j, k), readGrid(i, j, k)));
        }
      }
    }
  }
}

void poisson_xyz(Grid3D& Ex, Grid3D& Ey, Grid3D& Ez, Grid3D& rho, Greenfb& gf)
{
 int NX=rho.get_NX(),NY=rho.get_NX(), NZ=rho.get_NZ();
 double dx=rho.get_dx(), dy=rho.get_dy(), dz=rho.get_dz();
 
 poisson_xy(rho.get_ghostl2(),gf);
 poisson_xy(rho.get_ghostr2(),gf);
 poisson_xy(rho[0],gf); 
 
 for(int i=0; i<NZ; i++)
 {  
  if (i < NZ-1) poisson_xy(rho[i+1],gf);   
  for(int j=1;j<NX-1;j++)
   for(int l=1;l<NY-1;l++)
    {
     Ex[i](j,l)=-(rho[i](j+1,l)-rho[i](j-1,l))/(2.0*dx);
     Ey[i](j,l)=-(rho[i](j,l+1)-rho[i](j,l-1))/(2.0*dy);
     
     if (i==0) Ez[i](j,l)=-(rho[i+1](j,l)-rho.get_ghostl(j,l))/(2.0*dz);
     else if (i==NZ-1) Ez[i](j,l)=-(rho.get_ghostr(j,l)-rho[i-1](j,l))/(2.0*dz);
     else Ez[i](j,l)=-(rho[i+1](j,l)-rho[i-1](j,l))/(2.0*dz);
    }
 }
}

void MeshFileManager::writeShapePoints(const Grid3D &grid, const QString &path)
{
    QVector<Point3D> points = grid.getShapePoints();

    QFile file(getRootPath() + path);
    if (!file.open(QIODevice::WriteOnly)) {
        return;
    }

    QTextStream stream(&file);
    stream  << "MeshVersionFormatted 1" << endl
            << "Dimension" << endl
            << "3" << endl;

    stream  << "Vertices" << endl
            << points.size() << endl;

    for (int i = 0; i < points.size(); ++i) {
        stream << points[i].format()
               << endl;
    }

    stream << "End" << endl;
}

Rect Grid3DAction::getGridRect() const
{
    Grid3D *g = (Grid3D*)_gridNodeTarget->getGrid();
    return g->getGridRect();
}

Vertex3F Grid3DAction::getVertex(const Point& position) const
{
    Grid3D *g = (Grid3D*)_target->getGrid();
    return g->getVertex(position);
}

 static T evaluateFormula(Grid3D<T>& g, std::string exp) {
   expression = exp;
   g.generate(evalulateVoxelExpression);
 }

int main(int argc, char *argv[])
{
  // timing variables
  double start_total, stop_total, start_step, stop_step;
  #ifdef CPU_TIME
  double stop_init, init_min, init_max, init_avg;
  double start_bound, stop_bound, bound_min, bound_max, bound_avg;
  double start_hydro, stop_hydro, hydro_min, hydro_max, hydro_avg;
  double init, bound, hydro;
  init = bound = hydro = 0;
  #endif //CPU_TIME

  // start the total time
  start_total = get_time();

  /* Initialize MPI communication */
  #ifdef MPI_CHOLLA
  InitializeChollaMPI(&argc, &argv);
  #endif /*MPI_CHOLLA*/

  // declare Cfl coefficient and initial inverse timestep
  Real C_cfl = 0.4; // CFL coefficient 0 < C_cfl < 0.5 
  Real dti = 0; // inverse time step, 1.0 / dt

  // input parameter variables
  char *param_file;
  struct parameters P;
  int nfile = 0; // number of output files
  Real outtime = 0; // current output time

  // read in command line arguments
  if (argc != 2)
  {
    chprintf("usage: %s <parameter_file>\n", argv[0]);
    chexit(0);
  } else {
    param_file = argv[1];
  }

  // create the grid
  Grid3D G;

  // read in the parameters
  parse_params (param_file, &P);
  // and output to screen
  chprintf ("Parameter values:  nx = %d, ny = %d, nz = %d, tout = %f, init = %s, boundaries = %d %d %d %d %d %d\n", 
    P.nx, P.ny, P.nz, P.tout, P.init, P.xl_bcnd, P.xu_bcnd, P.yl_bcnd, P.yu_bcnd, P.zl_bcnd, P.zu_bcnd);
  chprintf ("Output directory:  %s\n", P.outdir);


  // initialize the grid
  G.Initialize(&P);
  chprintf("Local number of grid cells: %d %d %d %d\n", G.H.nx_real, G.H.ny_real, G.H.nz_real, G.H.n_cells);

  // Set initial conditions and calculate first dt
  chprintf("Setting initial conditions...\n");
  G.Set_Initial_Conditions(P, C_cfl);
  chprintf("Dimensions of each cell: dx = %f dy = %f dz = %f\n", G.H.dx, G.H.dy, G.H.dz);
  chprintf("Ratio of specific heats gamma = %f\n",gama);
  chprintf("Nstep = %d  Timestep = %f  Simulation time = %f\n", G.H.n_step, G.H.dt, G.H.t);
  // set main variables for Read_Grid inital conditions
  if (strcmp(P.init, "Read_Grid") == 0) {
    outtime += G.H.t;
    nfile = P.nfile;
    dti = 1.0 / G.H.dt;
  }

  // set boundary conditions (assign appropriate values to ghost cells)
  chprintf("Setting boundary conditions...\n");
  G.Set_Boundary_Conditions(P);
  chprintf("Boundary conditions set.\n");


  #ifdef OUTPUT
  // write the initial conditions to file
  chprintf("Writing initial conditions to file...\n");
  WriteData(G, P, nfile);
  // add one to the output file count
  nfile++;
  #endif //OUTPUT
  // increment the next output time
  outtime += P.outstep;

  #ifdef CPU_TIME
  stop_init = get_time();
  init = stop_init - start_total;
  #ifdef MPI_CHOLLA
  init_min = ReduceRealMin(init);
  init_max = ReduceRealMax(init);
  init_avg = ReduceRealAvg(init);
  chprintf("Init  min: %9.4f  max: %9.4f  avg: %9.4f\n", init_min, init_max, init_avg);
  #else
  printf("Init %9.4f\n", init);
  #endif //MPI_CHOLLA
  #endif //CPU_TIME


  // Evolve the grid, one timestep at a time
  chprintf("Starting caclulations.\n");
  while (G.H.t < P.tout)
//.........这里部分代码省略.........

//:
// Assumes given a Lookup table of integer squares.
// Also assumes the image \a im already has infinity in all non-zero points.
bool SaitoSquaredDistanceTransform::SDT_2D(Grid3D<GridElement>& grid, size_t sliceIndex, const std::vector<GridElement>& sq)
{
	const Tuple3ui& gridSize = grid.size();
	size_t r = gridSize.y;
	size_t c = gridSize.x;
	size_t voxelCount = r*c;

	GridElement* sliceData = grid.data() + sliceIndex * voxelCount;

	// 1st step: vertical row-wise EDT
	{
		if (!EDT_1D(sliceData, r, c))
		{
			return false;
		}
	}
	
	// 2nd step: horizontal scan
	{
		std::vector<GridElement> colData;
		try
		{
			colData.resize(r);
		}
		catch (const std::bad_alloc&)
		{
			//not enough memory
			return false;
		}

		//for each column
		for (size_t i = 0; i < c; ++i)
		{
			//fill buffer with column values
			{
				GridElement* pt = sliceData + i;
				for (size_t j = 0; j < r; ++j, pt += c)
					colData[j] = *pt;
			}

			//forward scan
			GridElement* pt = sliceData + i + c;
			{
				GridElement a = 0;
				GridElement buffer = colData[0];

				for (size_t j = 1; j < r; ++j, pt += c)
				{
					size_t rowIndex = j;
					if (a != 0)
						--a;
					if (colData[rowIndex] > buffer + 1)
					{
						GridElement b = (colData[rowIndex] - buffer - 1) / 2;
						if (rowIndex + b + 1 > r)
							b = static_cast<GridElement>(r - 1 - rowIndex);

						GridElement* npt = pt + a*c;
						for (GridElement l = a; l <= b; ++l)
						{
							GridElement m = buffer + sq[l + 1];
							if (colData[rowIndex + l] <= m)
							{
								//proceed to next column
								break;
							}
							if (m < *npt)
								*npt = m;
							npt += c;
						}
						a = b;
					}
					else
					{
						a = 0;
					}
					buffer = colData[rowIndex];
				}
			}

			//backward scan
			pt -= 2 * c;
			{
				GridElement a = 0;
				GridElement buffer = colData[r - 1];

				for (size_t j = 1; j < r; ++j, pt -= c)
				{
					size_t rowIndex = r - j - 1;
					if (a != 0)
						--a;
					if (colData[rowIndex] > buffer + 1)
					{
						GridElement b = (colData[rowIndex] - buffer - 1) / 2;
						if (rowIndex < b)
							b = static_cast<GridElement>(rowIndex);

//.........这里部分代码省略.........

bool SaitoSquaredDistanceTransform::SDT_3D(Grid3D<GridElement>& grid, GenericProgressCallback* progressCb/*=0*/)
{
	const Tuple3ui& gridSize = grid.size();
	size_t r = gridSize.y;
	size_t c = gridSize.x;
	size_t p = gridSize.z;
	size_t voxelCount = r*c*p;

	size_t diag = static_cast<size_t>(ceil(sqrt(static_cast<double>(r*r + c*c + p*p))) - 1);
	size_t nsqr = 2 * (diag + 1);

	std::vector<GridElement> sq;
	try
	{
		sq.resize(nsqr);
	}
	catch (const std::bad_alloc&)
	{
		//not enough memory
		return false;
	}

	for (size_t i = 0; i < nsqr; ++i)
	{
		sq[i] = static_cast<GridElement>(i*i);
	}

	const GridElement maxDistance = std::numeric_limits<GridElement>::max() - static_cast<GridElement>(r*r - c*c - p*p) - 1;

	NormalizedProgress normProgress(progressCb, static_cast<unsigned>(p + r));
	if (progressCb)
	{
		progressCb->setMethodTitle("Saito Distance Transform");
		char buffer[256];
		sprintf(buffer, "Box: [%u x %u x %u]", gridSize.x, gridSize.y, gridSize.z);
		progressCb->setInfo(buffer);
		progressCb->reset();
		progressCb->start();
	}

	GridElement* data = grid.data();
	{
		for (size_t i = 0; i < voxelCount; ++i)
		{
			//DGM: warning we must invert the input image here!
			if (data[i] == 0)
				data[i] = maxDistance;
			else
				data[i] = 0;
		}
	}

	// 2D EDT for each slice
	for (size_t k = 0; k < p; ++k)
	{
		if (!SDT_2D(grid, k, sq))
		{
			return false;
		}

		if (progressCb && !normProgress.oneStep())
		{
			//process cancelled by user
			return false;
		}
	}

	// Now, for each pixel, compute final distance by searching along Z direction
	size_t rc = r*c;
	std::vector<GridElement> colData;
	try
	{
		colData.resize(p);
	}
	catch (const std::bad_alloc&)
	{
		//not enough memory
		return false;
	}

	for (size_t j = 0; j < r; ++j, data += c)
	{
		for (size_t i = 0; i < c; ++i)
		{
			GridElement* pt = data + i;

			for (size_t k = 0; k < p; ++k, pt += rc)
				colData[k] = *pt;

			pt = data + i + rc;
			GridElement a = 0;
			GridElement buffer = colData[0];

			for (size_t k = 1; k < p; ++k, pt += rc)
			{
				if (a != 0)
					--a;
				if (colData[k] > buffer + 1)
				{
					GridElement b = (colData[k] - buffer - 1) / 2;
//.........这里部分代码省略.........

void Pic::gatherXYZ(double pic_charge, Grid3D& target){
  long j;
  target.reset();
  for(j=0; j<pics.size(); ++j)
    target.Pic2Grid(pic_charge, pics[j].x, pics[j].y, pics[j].z);
}