C++ Index::X方法代码示例

GlobalGridPartitioning::GlobalGridPartitioning(Index pos, Index procs, Index size, int points_min)
{
  Index i;
  Index index_beg, index_end;
  Index local_size;

  const Index active_procs = size / points_min;
  const Index size_per_proc = size / active_procs;
  const Index remainder = size % active_procs;

  for (i.X() = 0; i.X() < active_procs.X(); ++i.X())
    for (i.Y() = 0; i.Y() < active_procs.Y(); ++i.Y())
      for (i.Z() = 0; i.Z() < active_procs.Z(); ++i.Z()) {

	local_size = size_per_proc;

	for (int j=0; j<3; ++j)
	  if (pos[j] < remainder[j])
	    ++(local_size[j]);

	index_beg = i * local_size;

	for (int j=0; j<3; ++j)
	  if (i[j] >= remainder[j])
	    index_beg[j] += remainder[j];

	index_end = index_beg + local_size;

	begin.insert(std::make_pair(i, index_beg));
	end.insert(std::make_pair(i, index_end));

      }

}

vmg_float Helper::InterpolateTrilinear(const Vector& point, const Grid& grid)
{
  vmg_float interpolate_vals[4], grid_vals[8];

  const Index index_global = (point - grid.Extent().Begin()) / grid.Extent().MeshWidth();
  const Index index_local = index_global - grid.Global().LocalBegin() + grid.Local().Begin();
  const Vector coord = (point - grid.Extent().Begin() - index_global * grid.Extent().MeshWidth()) / grid.Extent().MeshWidth();

  grid_vals[0] = grid.GetVal(index_local.X()  , index_local.Y()  , index_local.Z()  );
  grid_vals[1] = grid.GetVal(index_local.X()+1, index_local.Y()  , index_local.Z()  );
  grid_vals[2] = grid.GetVal(index_local.X()  , index_local.Y()+1, index_local.Z()  );
  grid_vals[3] = grid.GetVal(index_local.X()+1, index_local.Y()+1, index_local.Z()  );
  grid_vals[4] = grid.GetVal(index_local.X()  , index_local.Y()  , index_local.Z()+1);
  grid_vals[5] = grid.GetVal(index_local.X()+1, index_local.Y()  , index_local.Z()+1);
  grid_vals[6] = grid.GetVal(index_local.X()  , index_local.Y()+1, index_local.Z()+1);
  grid_vals[7] = grid.GetVal(index_local.X()+1, index_local.Y()+1, index_local.Z()+1);

  for (int i=0; i<4; ++i)
    interpolate_vals[i] = (1.0 - coord.X()) * grid_vals[2*i] + coord.X() * grid_vals[2*i+1];

  for (int i=0; i<2; ++i)
    interpolate_vals[i] = (1.0 - coord.Y()) * interpolate_vals[2*i] + coord.Y() * interpolate_vals[2*i+1];

  return (1.0 - coord.Z()) * interpolate_vals[0] + coord.Z() * interpolate_vals[1];
}

static inline void ComputePartial(Grid& sol, Grid& rhs, const Stencil& mat,
				  const Index& begin, const Index& end,
				  const vmg_float& prefactor, const vmg_float& diag_inv,
				  const int& off)
{
  int i,j,k;
  vmg_float temp;
  Stencil::iterator iter;

  for (i=begin.X(); i<end.X(); ++i)
    for (j=begin.Y(); j<end.Y(); ++j) {

      int z_begin = begin.Z() + (i + j + begin.Z() + off) % 2;

#ifdef DEBUG
      int off_sum = MG::GetComm()->LevelSum(rhs, z_begin - begin.Z());
      assert(z_begin - begin.Z() == 0 || z_begin - begin.Z() == 1);
      assert(off_sum == 0 || off_sum == MG::GetComm()->Size(rhs));
#endif /* DEBUG */

      for (k=z_begin; k<end.Z(); k+=2) {

	temp = prefactor * rhs.GetVal(i,j,k);

	for (iter=mat.begin(); iter!=mat.end(); ++iter)
	  temp -= iter->Val() * sol.GetVal(i+iter->Disp().X(),
					   j+iter->Disp().Y(),
					   k+iter->Disp().Z());

	sol(i,j,k) = temp * diag_inv;

      }
    }
}

static inline void ComputePartial(Grid& sol, Grid& rhs,
				  const Index& begin, const Index& end,
				  const vmg_float& prefactor, const int& off)
{
  const vmg_float fac = 1.0 / 6.0;

  for (int i=begin.X(); i<end.X(); ++i)
    for (int j=begin.Y(); j<end.Y(); ++j)
      for (int k=begin.Z() + (i + j + begin.Z() + off) % 2; k<end.Z(); k+=2)
	sol(i,j,k) = prefactor * rhs.GetVal(i,j,k) + fac * (sol.GetVal(i-1,j  ,k  ) +
							    sol.GetVal(i+1,j  ,k  ) +
							    sol.GetVal(i  ,j-1,k  ) +
							    sol.GetVal(i  ,j+1,k  ) +
							    sol.GetVal(i  ,j  ,k-1) +
							    sol.GetVal(i  ,j  ,k+1));
}

void CommSerial::PrintGrid(Grid& grid, const char* information)
{
  Index i;
  std::stringstream out;
  std::ofstream out_file;

  OpenFileAndPrintHeader(out_file, grid, information);

  for (i.Z()=grid.Local().Begin().Z(); i.Z()<grid.Local().End().Z(); ++i.Z())
    for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
      for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
	out << std::scientific << grid.GetVal(i) << std::endl;

  out_file << out.str();
  out_file.close();
}

void InterfaceParticles::ExportSolution(Grid& grid)
{
  Index i;

#ifdef OUTPUT_DEBUG
  vmg_float e = 0.0;
  vmg_float e_long = 0.0;
  vmg_float e_self = 0.0;
  vmg_float e_short_peak = 0.0;
  vmg_float e_short_spline = 0.0;
#endif

  Factory& factory = MG::GetFactory();
  Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());

  /*
   * Get parameters and arrays
   */
  const vmg_int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
  const vmg_int& interpolation_degree = factory.GetObjectStorageVal<int>("PARTICLE_INTERPOLATION_DEGREE");

  Particle::Interpolation ip(interpolation_degree);

  const vmg_float r_cut = near_field_cells * grid.Extent().MeshWidth().Max();

  /*
   * Copy potential values to a grid with sufficiently large halo size.
   * This may be optimized in future.
   * The parameters of this grid have been set in the import step.
   */
  Grid& particle_grid = comm.GetParticleGrid();

  for (i.X()=0; i.X()<grid.Local().Size().X(); ++i.X())
    for (i.Y()=0; i.Y()<grid.Local().Size().Y(); ++i.Y())
      for (i.Z()=0; i.Z()<grid.Local().Size().Z(); ++i.Z())
        particle_grid(i + particle_grid.Local().Begin()) = grid.GetVal(i + grid.Local().Begin());

  comm.CommToGhosts(particle_grid);

  /*
   * Compute potentials
   */
  Particle::LinkedCellList lc(particles, near_field_cells, grid);
  Particle::LinkedCellList::iterator p1, p2;
  Grid::iterator iter;

  comm.CommLCListToGhosts(lc);

  for (int i=lc.Local().Begin().X(); i<lc.Local().End().X(); ++i)
    for (int j=lc.Local().Begin().Y(); j<lc.Local().End().Y(); ++j)
      for (int k=lc.Local().Begin().Z(); k<lc.Local().End().Z(); ++k) {

	if (lc(i,j,k).size() > 0)
	  ip.ComputeCoefficients(particle_grid, Index(i,j,k) - lc.Local().Begin() + particle_grid.Local().Begin());

	for (p1=lc(i,j,k).begin(); p1!=lc(i,j,k).end(); ++p1) {

	  // Interpolate long-range part of potential and electric field
	  ip.Evaluate(**p1);

	  // Subtract self-induced potential
	  (*p1)->Pot() -= (*p1)->Charge() * spl.GetAntiDerivativeAtZero();

#ifdef OUTPUT_DEBUG
	  e_long += 0.5 * (*p1)->Charge() * ip.EvaluatePotentialLR(**p1);
	  e_self += 0.5 * (*p1)->Charge() * (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
#endif

	  for (int dx=-1*near_field_cells; dx<=near_field_cells; ++dx)
	    for (int dy=-1*near_field_cells; dy<=near_field_cells; ++dy)
	      for (int dz=-1*near_field_cells; dz<=near_field_cells; ++dz) {

		for (p2=lc(i+dx,j+dy,k+dz).begin(); p2!=lc(i+dx,j+dy,k+dz).end(); ++p2)

		  if (*p1 != *p2) {

		    const Vector dir = (*p1)->Pos() - (*p2)->Pos();
		    const vmg_float length = dir.Length();

		    if (length < r_cut) {

		      (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
		      (*p1)->Field() += (*p2)->Charge() * dir * spl.EvaluateField(length);

#ifdef OUTPUT_DEBUG
		      e_short_peak += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
		      e_short_spline += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
#endif
		    }
		  }
	      }
	}
      }

  /* Remove average force term */
  Vector average_force = 0.0;
  for (std::list<Particle::Particle>::const_iterator iter=particles.begin(); iter!=particles.end(); ++iter)
    average_force += iter->Charge() * iter->Field();
  const vmg_int& npl = MG::GetFactory().GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
  const vmg_int num_particles_global = comm.GlobalSum(npl);
//.........这里部分代码省略.........

inline int Grid::GlobalLinearIndex(const Index& index) const
{
  return GlobalLinearIndex(index.X(), index.Y(), index.Z());
}

inline int IsGrid<T>::GlobalLinearIndex(const Index& index) const
{
  return index.Z() + global.GlobalSize().Z() * (index.Y() + global.GlobalSize().Y() * index.X());
}

void Particle::CommMPI::CommLCListToGhosts(LinkedCellList& lc)
{
  VMG::MPI::DatatypesLocal types(lc, comm_global, false);
  std::vector<int> send_size(types.NB().size());
  vmg_int recv_size;
  std::list<Particle*>::iterator iter;
  Index ind;
  Vector offset;

  const Vector halo_length = lc.Local().HaloSize1() * lc.Extent().MeshWidth();

  lc.ClearHalo();

  for (unsigned int i=0; i<types.NB().size(); ++i)
    if (types.NB()[i].Feasible()) {

      for (int j=0; j<3; ++j)
        if ((types.Offset()[i][j] < 0 && lc.Global().LocalBegin()[j] == 0) ||
            (types.Offset()[i][j] > 0 && lc.Global().LocalEnd()[j] == lc.Global().GlobalSize()[j]))
          offset[j] = -1.0 * types.Offset()[i][j] * lc.Extent().Size()[j];
        else
          offset[j] = 0.0;

      for (ind.X() = types.NB()[i].Starts().X(); ind.X() < types.NB()[i].Starts().X()+types.NB()[i].Subsizes().X(); ++ind.X())
	for (ind.Y() = types.NB()[i].Starts().Y(); ind.Y() < types.NB()[i].Starts().Y()+types.NB()[i].Subsizes().Y(); ++ind.Y())
	  for (ind.Z() = types.NB()[i].Starts().Z(); ind.Z() < types.NB()[i].Starts().Z()+types.NB()[i].Subsizes().Z(); ++ind.Z())
	    for (iter=lc(ind).begin(); iter!=lc(ind).end(); ++iter) {

              for (int j=0; j<3; ++j)
                types.NB()[i].Buffer().push_back((*iter)->Pos()[j] + offset[j]);
              types.NB()[i].Buffer().push_back((*iter)->Charge());

              assert(lc.Extent().Begin().IsComponentwiseLessOrEqual((*iter)->Pos()));
              assert(lc.Extent().End().IsComponentwiseGreaterOrEqual((*iter)->Pos()));
              assert(lc.Extent().Begin().IsComponentwiseLessOrEqual((*iter)->Pos() + offset + halo_length));
              assert(lc.Extent().End().IsComponentwiseGreaterOrEqual((*iter)->Pos() + offset - halo_length));
	    }

      send_size[i] = types.NB()[i].Buffer().size();
      MPI_Isend(&send_size[i], 1, MPI_INT, types.NB()[i].Rank(), 2048+types.NB()[i].TagSend(), comm_global, &Request());

      if (send_size[i] > 0)
	MPI_Isend(&types.NB()[i].Buffer().front(), send_size[i], MPI_DOUBLE,
		  types.NB()[i].Rank(), 4096+types.NB()[i].TagSend(),
		  comm_global, &Request());
    }

  for (unsigned int i=0; i<types.Halo().size(); ++i)
    if (types.Halo()[i].Feasible()) {
      MPI_Recv(&recv_size, 1, MPI_INT, types.Halo()[i].Rank(), 2048+types.Halo()[i].TagReceive(), comm_global, MPI_STATUS_IGNORE);
      if (recv_size > 0) {
	types.Halo()[i].Buffer().resize(recv_size);
	MPI_Irecv(&types.Halo()[i].Buffer().front(), recv_size, MPI_DOUBLE,
		  types.Halo()[i].Rank(), 4096+types.Halo()[i].TagReceive(),
		  comm_global, &Request());
      }
    }

  WaitAll();

  for (unsigned int i=0; i<types.Halo().size(); ++i)
    for (unsigned int j=0; j<types.Halo()[i].Buffer().size(); j+=4)
      lc.AddParticleToHalo(&types.Halo()[i].Buffer()[j], types.Halo()[i].Buffer()[j+3]);
}