C++ ParallelComm类代码示例

void test_read_parallel(int num_verts)
{
  Core moab;
  Interface& mb = moab;
  EntityHandle file_set;
  ErrorCode rval;
  rval = mb.create_meshset(MESHSET_SET, file_set);
  CHECK_ERR(rval);

  std::string opt = std::string("PARALLEL=READ_PART;PARTITION=;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS") +
      partition_method;
  rval = mb.load_file(example, &file_set, opt.c_str());
  CHECK_ERR(rval);

  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);

  rval = pcomm->check_all_shared_handles();
  CHECK_ERR(rval);

    // get the total # owned verts
  Range verts;
  rval = mb.get_entities_by_type(0, MBVERTEX, verts);
  CHECK_ERR(rval);
  rval = pcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT);
  CHECK_ERR(rval);
  int my_num = verts.size(), total_verts;
  MPI_Reduce(&my_num, &total_verts, 1, MPI_INTEGER, MPI_SUM, 0, pcomm->proc_config().proc_comm());
  
  if (0 == pcomm->proc_config().proc_rank()) CHECK_EQUAL(total_verts, num_verts);
}

void test_read_conn()
{
  Core moab;
  Interface& mb = moab;

  std::string opts;
  get_options(opts);

  ErrorCode rval = mb.load_file(conn_fname, NULL, opts.c_str());
  CHECK_ERR(rval);

#ifdef MOAB_HAVE_MPI
  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
#else
  int procs = 1;
#endif

  // Make check runs this test on one processor
  if (1 == procs) {
    // Get vertices
    Range verts;
    rval = mb.get_entities_by_type(0, MBVERTEX, verts);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)3458, verts.size());

    // Get cells
    Range cells;
    rval = mb.get_entities_by_type(0, MBQUAD, cells);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)3456, cells.size());
  }
}

void intersection_at_level(iMesh_Instance instance,
    iBase_EntitySetHandle fine_set, iBase_EntitySetHandle lagr_set,
    iBase_EntitySetHandle intx_set, double * dep_coords, double radius2,
    int * ierr) {
  *ierr = 1;
  Interface * mb = MOABI;
  //MPI_Comm mpicomm = MPI_Comm_f2c(comm);
  // instantiate parallel comm now or not?

  EntityHandle lagrMeshSet = (EntityHandle) lagr_set;

  ParallelComm *pcomm = ParallelComm::get_pcomm(mb, 0);
  if (NULL == pcomm)
    return; // error is 1

  // set the departure tag on the fine mesh vertices
  ErrorCode rval = set_departure_points_position(mb, lagrMeshSet, dep_coords,
      radius2);
  ERRORV(rval, "can't set departure tag");
  if (debug) {
    std::stringstream fff;
    fff << "lagr0" << pcomm->proc_config().proc_rank() << ".vtk";
    rval = mb->write_mesh(fff.str().c_str(), &lagrMeshSet, 1);
    ERRORV(rval, "can't write covering set ");
  }

  // it should be done earlier
  pworker->SetRadius(radius);

  EntityHandle covering_set;
  rval = pworker->create_departure_mesh_3rd_alg(lagrMeshSet, covering_set);
  ERRORV(rval, "can't compute covering set ");

  if (debug) {
    std::stringstream fff;
    fff << "cover" << pcomm->proc_config().proc_rank() << ".vtk";
    rval = mb->write_mesh(fff.str().c_str(), &covering_set, 1);

    ERRORV(rval, "can't write covering set ");
  }
  EntityHandle intxSet = (EntityHandle) intx_set;
  rval = pworker->intersect_meshes(covering_set, (EntityHandle) fine_set,
      intxSet);
  ERRORV(rval, "can't intersect ");

  if (debug) {
    std::stringstream fff;
    fff << "intx0" << pcomm->proc_config().proc_rank() << ".vtk";
    rval = mb->write_mesh(fff.str().c_str(), &intxSet, 1);
    ERRORV(rval, "can't write covering set ");
  }

  return;
}

void test_read_onevar() 
{
  Core moab;
  Interface& mb = moab;

  std::string opts;
  get_options(opts);

  // Read mesh and read vertex variable T at all timesteps
  opts += std::string(";VARIABLE=T");
  ErrorCode rval = mb.load_file(example, NULL, opts.c_str());
  CHECK_ERR(rval);

#ifdef MOAB_HAVE_MPI
  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
#else
  int procs = 1;
#endif

  // Make check runs this test on one processor
  if (1 == procs) {
    // Check for proper tags
    Tag Ttag0, Ttag1;
    rval = mb.tag_get_handle("T0", levels, MB_TYPE_DOUBLE, Ttag0);
    CHECK_ERR(rval);
    rval = mb.tag_get_handle("T1", levels, MB_TYPE_DOUBLE, Ttag1);
    CHECK_ERR(rval);

    // Get vertices
    Range verts;
    rval = mb.get_entities_by_type(0, MBVERTEX, verts);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)3458, verts.size());

    // Get all values of tag T0
    int count;
    void* Tbuf;
    rval = mb.tag_iterate(Ttag0, verts.begin(), verts.end(), count, Tbuf);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)count, verts.size());

    const double eps = 0.0001;
    double* data = (double*) Tbuf;

    // Check first level values on 4 strategically selected vertices
    CHECK_REAL_EQUAL(233.1136, data[0 * levels], eps); // First vert
    CHECK_REAL_EQUAL(236.1505, data[1728 * levels], eps); // Median vert
    CHECK_REAL_EQUAL(235.7722, data[1729 * levels], eps); // Median vert
    CHECK_REAL_EQUAL(234.0416, data[3457 * levels], eps); // Last vert
  }
}

ErrorCode test_read(const char *filename, const char *option) 
{
  Core mb_instance;
  Interface& moab = mb_instance;
  ErrorCode rval;

  rval = moab.load_file( filename, 0, option);
  CHKERR(rval);

  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab, 0);

  rval = pcomm->check_all_shared_handles();
  CHKERR(rval);

  return MB_SUCCESS;
}

/**
  Perform full parallel sweep algorithm on subset of subdomains.
*/
void SweepSubdomains (std::vector<int> subdomain_list, Grid_Data *grid_data, bool block_jacobi)
{
  // Create a new sweep communicator object
  ParallelComm *comm = NULL;
  if(block_jacobi){
    comm = new BlockJacobiComm(grid_data);
  }
  else {
    comm = new SweepComm(grid_data);
  }

  // Add all subdomains in our list
  for(int i = 0;i < subdomain_list.size();++ i){
    int sdom_id = subdomain_list[i];
    comm->addSubdomain(sdom_id, grid_data->subdomains[sdom_id]);
  }

  /* Loop until we have finished all of our work */
  while(comm->workRemaining()){

    // Get a list of subdomains that have met dependencies
    std::vector<int> sdom_ready = comm->readySubdomains();
    int backlog = sdom_ready.size();

    // Run top of list
    if(backlog > 0){
      int sdom_id = sdom_ready[0];
      Subdomain &sdom = grid_data->subdomains[sdom_id];
      // Clear boundary conditions
      for(int dim = 0;dim < 3;++ dim){
        if(sdom.upwind[dim].subdomain_id == -1){
          sdom.plane_data[dim]->clear(0.0);
        }
      }
      {
        BLOCK_TIMER(grid_data->timing, Sweep_Kernel);
        // Perform subdomain sweep
        grid_data->kernel->sweep(&sdom);
      }

      // Mark as complete (and do any communication)
      comm->markComplete(sdom_id);
    }
  }

  delete comm;
}

void multiple_loads_of_same_file()
{
  Core moab;
  Interface& mb = moab;

  // Need a file set for nomesh to work right
  EntityHandle file_set;
  ErrorCode rval;
  rval = mb.create_meshset(MESHSET_SET, file_set);
  CHECK_ERR(rval);

  // Read first only header information, no mesh, no variable
  read_options = "PARALLEL=READ_PART;PARTITION;NOMESH;VARIABLE=;PARTITION_METHOD=TRIVIAL";

  rval = mb.load_file(example, &file_set, read_options.c_str());
  CHECK_ERR(rval);

  // Create mesh, no variable
  read_options = "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION_METHOD=TRIVIAL;VARIABLE=";

  rval = mb.load_file(example, &file_set, read_options.c_str());
  CHECK_ERR(rval);

  // Read variable vorticity at timestep 0, no mesh
  read_options = "PARALLEL=READ_PART;PARTITION;PARTITION_METHOD=TRIVIAL;NOMESH;VARIABLE=vorticity;TIMESTEP=0";

  rval = mb.load_file(example, &file_set, read_options.c_str());
  CHECK_ERR(rval);

  Range local_verts;
  rval = mb.get_entities_by_type(file_set, MBVERTEX, local_verts);
  CHECK_ERR(rval);

  Range local_edges;
  rval = mb.get_entities_by_type(file_set, MBEDGE, local_edges);
  CHECK_ERR(rval);

  Range local_cells;
  rval = mb.get_entities_by_type(file_set, MBPOLYGON, local_cells);
  CHECK_ERR(rval);
  // No mixed elements
  CHECK_EQUAL((size_t)1, local_cells.psize());

  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
  int rank = pcomm->proc_config().proc_rank();

  // Make check runs this test on two processors
  if (2 == procs) {
    CHECK_EQUAL((size_t)321, local_cells.size());

    // Check tag for cell variable vorticity at timestep 0
    Tag vorticity_tag0;
    rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0);
    CHECK_ERR(rval);

    // Get vorticity0 tag values on 3 local cells
    double vorticity0_val[3 * layers];
    EntityHandle cell_ents[] = {local_cells[0], local_cells[160], local_cells[320]};
    rval = mb.tag_get_data(vorticity_tag0, cell_ents, 3, vorticity0_val);
    CHECK_ERR(rval);

    if (0 == rank) {
      CHECK_EQUAL((size_t)687, local_verts.size());
      CHECK_EQUAL((size_t)1007, local_edges.size());

      // Layer 0
      CHECK_REAL_EQUAL(3.629994, vorticity0_val[0 * layers], eps);
      CHECK_REAL_EQUAL(-1.708188, vorticity0_val[1 * layers], eps);
      CHECK_REAL_EQUAL(0.131688, vorticity0_val[2 * layers], eps);
      // Layer 1
      CHECK_REAL_EQUAL(3.629944, vorticity0_val[0 * layers + 1], eps);
      CHECK_REAL_EQUAL(-1.708164, vorticity0_val[1 * layers + 1], eps);
      CHECK_REAL_EQUAL(0.131686, vorticity0_val[2 * layers + 1], eps);
    }
    else if (1 == rank) {
      CHECK_EQUAL((size_t)688, local_verts.size());
      CHECK_EQUAL((size_t)1008, local_edges.size());

      // Layer 0
      CHECK_REAL_EQUAL(-0.554888, vorticity0_val[0 * layers], eps);
      CHECK_REAL_EQUAL(2.434397, vorticity0_val[1 * layers], eps);
      CHECK_REAL_EQUAL(-0.554888, vorticity0_val[2 * layers], eps);
      // Layer 1
      CHECK_REAL_EQUAL(-0.554881, vorticity0_val[0 * layers + 1], eps);
      CHECK_REAL_EQUAL(2.434363, vorticity0_val[1 * layers + 1], eps);
      CHECK_REAL_EQUAL(-0.554881, vorticity0_val[2 * layers + 1], eps);
    }
  }
}

void read_mesh_parallel(bool rcbzoltan)
{
  Core moab;
  Interface& mb = moab;

  read_options = "PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=";
  if (rcbzoltan)
    read_options = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=";

  ErrorCode rval = mb.load_file(example, NULL, read_options.c_str());
  CHECK_ERR(rval);

  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
  int rank = pcomm->proc_config().proc_rank();

  rval = pcomm->check_all_shared_handles();
  CHECK_ERR(rval);

  // Get local vertices
  Range local_verts;
  rval = mb.get_entities_by_type(0, MBVERTEX, local_verts);
  CHECK_ERR(rval);

  int verts_num = local_verts.size();
  if (2 == procs) {
    if (rcbzoltan) {
      if (0 == rank)
        CHECK_EQUAL(684, verts_num);
      else if (1 == rank)
        CHECK_EQUAL(691, verts_num); // Not owned vertices included
    }
    else {
      if (0 == rank)
        CHECK_EQUAL(687, verts_num);
      else if (1 == rank)
        CHECK_EQUAL(688, verts_num); // Not owned vertices included
    }
  }

  rval = pcomm->filter_pstatus(local_verts, PSTATUS_NOT_OWNED, PSTATUS_NOT);
  CHECK_ERR(rval);

  verts_num = local_verts.size();
  if (2 == procs) {
    if (rcbzoltan) {
      if (0 == rank)
        CHECK_EQUAL(684, verts_num);
      else if (1 == rank)
        CHECK_EQUAL(596, verts_num); // Not owned vertices excluded
    }
    else {
      if (0 == rank)
        CHECK_EQUAL(687, verts_num);
      else if (1 == rank)
        CHECK_EQUAL(593, verts_num); // Not owned vertices excluded
    }
  }

  // Get local edges
  Range local_edges;
  rval = mb.get_entities_by_type(0, MBEDGE, local_edges);
  CHECK_ERR(rval);

  int edges_num = local_edges.size();
  if (2 == procs) {
    if (rcbzoltan) {
      if (0 == rank)
        CHECK_EQUAL(1002, edges_num);
      else if (1 == rank)
        CHECK_EQUAL(1013, edges_num); // Not owned edges included
    }
    else {
      if (0 == rank)
        CHECK_EQUAL(1007, edges_num);
      else if (1 == rank)
        CHECK_EQUAL(1008, edges_num); // Not owned edges included
    }
  }

  rval = pcomm->filter_pstatus(local_edges, PSTATUS_NOT_OWNED, PSTATUS_NOT);
  CHECK_ERR(rval);

  edges_num = local_edges.size();
  if (2 == procs) {
    if (rcbzoltan) {
      if (0 == rank)
        CHECK_EQUAL(1002, edges_num);
      else if (1 == rank)
        CHECK_EQUAL(918, edges_num); // Not owned edges excluded
    }
    else {
      if (0 == rank)
        CHECK_EQUAL(1007, edges_num);
      else if (1 == rank)
        CHECK_EQUAL(913, edges_num); // Not owned edges excluded
    }
  }

  // Get local cells
//.........这里部分代码省略.........

void gather_one_cell_var(int gather_set_rank)
{
  Core moab;
  Interface& mb = moab;

  EntityHandle file_set;
  ErrorCode rval = mb.create_meshset(MESHSET_SET, file_set);
  CHECK_ERR(rval);

  read_options = "PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS";
  std::ostringstream gather_set_option;
  gather_set_option << ";GATHER_SET=" << gather_set_rank;
  read_options += gather_set_option.str();

  rval = mb.load_file(example, &file_set, read_options.c_str());
  CHECK_ERR(rval);

  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
  int rank = pcomm->proc_config().proc_rank();

  // Make sure gather_set_rank is valid
  if (gather_set_rank < 0 || gather_set_rank >= procs)
    return;

  Range cells, cells_owned;
  rval = mb.get_entities_by_type(file_set, MBPOLYGON, cells);
  CHECK_ERR(rval);

  // Get local owned cells
  rval = pcomm->filter_pstatus(cells, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &cells_owned);
  CHECK_ERR(rval);

  EntityHandle gather_set = 0;
  if (gather_set_rank == rank) {
    // Get gather set
    ReadUtilIface* readUtilIface;
    mb.query_interface(readUtilIface);
    rval = readUtilIface->get_gather_set(gather_set);
    CHECK_ERR(rval);
    assert(gather_set != 0);
  }

  Tag vorticity_tag0, gid_tag;
  rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0, MB_TAG_DENSE);
  CHECK_ERR(rval);

  rval = mb.tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid_tag, MB_TAG_DENSE);
  CHECK_ERR(rval);

  pcomm->gather_data(cells_owned, vorticity_tag0, gid_tag, gather_set, gather_set_rank);

  if (gather_set_rank == rank) {
    // Get gather set cells
    Range gather_set_cells;
    rval = mb.get_entities_by_type(gather_set, MBPOLYGON, gather_set_cells);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)642, gather_set_cells.size());
    CHECK_EQUAL((size_t)1, gather_set_cells.psize());

    // Check vorticity0 tag values on 4 gather set cells: first cell, two median cells, and last cell
    EntityHandle cell_ents[] = {gather_set_cells[0], gather_set_cells[320],
                                gather_set_cells[321], gather_set_cells[641]};
    double vorticity0_val[4 * layers];
    rval = mb.tag_get_data(vorticity_tag0, &cell_ents[0], 4, vorticity0_val);
    CHECK_ERR(rval);

    // Only check first two layers
    // Layer 0
    CHECK_REAL_EQUAL(3.629994, vorticity0_val[0 * layers], eps);
    CHECK_REAL_EQUAL(0.131688, vorticity0_val[1 * layers], eps);
    CHECK_REAL_EQUAL(-0.554888, vorticity0_val[2 * layers], eps);
    CHECK_REAL_EQUAL(-0.554888, vorticity0_val[3 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(3.629944, vorticity0_val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(0.131686, vorticity0_val[1 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.554881, vorticity0_val[2 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.554881, vorticity0_val[3 * layers + 1], eps);
  }
}

int main(int argc, char **argv)
{

  MPI_Init(&argc, &argv);
  LONG_DESC << "This program simulates a transport problem on a sphere"
        " according to a benchmark from a Nair & Lauritzen paper.\n"
        << "It starts with a partitioned mesh on a sphere, add a tracer, and steps through.\n" <<
        "The flow reverses after half time, and it should return to original configuration, if the integration was exact. ";
  ProgOptions opts(LONG_DESC.str(), BRIEF_DESC);

  // read a homme file, partitioned in 16 so far
  std::string fileN= TestDir + "/HN16.h5m";
  const char *filename_mesh1 = fileN.c_str();

  opts.addOpt<double>("gtolerance,g",
      "geometric absolute tolerance (used for point concidence on the sphere)", &gtol);

  std::string input_file;
  opts.addOpt<std::string>("input_file,i", "input mesh file, partitioned",
      &input_file);
  std::string extra_read_opts;
  opts.addOpt<std::string>("extra_read_options,O", "extra read options ",
        &extra_read_opts);
  //int field_type;
  opts.addOpt<int>("field_type,f",
        "field type--  1: quasi-smooth; 2: smooth; 3: slotted cylinders (non-smooth)", &field_type);

  opts.addOpt<int>("num_steps,n",
          "number of  steps ", &numSteps);

  //bool reorder = false;
  opts.addOpt<void>("write_debug_files,w", "write debugging files during simulation ",
        &writeFiles);

  opts.addOpt<void>("write_velocity_files,v", "Reorder mesh to group entities by partition",
     &velocity);

  opts.addOpt<void>("write_result_in_parallel,p", "write tracer result files",
     &parallelWrite);

  opts.parseCommandLine(argc, argv);

  if (!input_file.empty())
    filename_mesh1=input_file.c_str();

  // read in parallel, in the "euler_set", the initial mesh
  std::string optsRead = std::string("PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION")+
            std::string(";PARALLEL_RESOLVE_SHARED_ENTS")+extra_read_opts;
  Core moab;
  Interface & mb = moab;
  EntityHandle euler_set;
  ErrorCode rval;
  rval = mb.create_meshset(MESHSET_SET, euler_set);
  CHECK_ERR(rval);

  rval = mb.load_file(filename_mesh1, &euler_set, optsRead.c_str());

  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  CHECK_ERR(rval);

  rval = pcomm->check_all_shared_handles();
  CHECK_ERR(rval);

  int rank = pcomm->proc_config().proc_rank();

  if (0==rank)
  {
    std::cout << " case 1: use -gtol " << gtol <<
        " -R " << radius << " -input " << filename_mesh1 <<  " -f " << field_type <<
        " numSteps: " << numSteps << "\n";
    std::cout<<" write debug results: " << (writeFiles ? "yes" : "no") << "\n";
    std::cout<< " write tracer in parallel: " << ( parallelWrite ? "yes" : "no") << "\n";
    std::cout <<" output velocity: " << (velocity? "yes" : "no") << "\n";
  }

  // tagTracer is the value at nodes
  Tag tagTracer = 0;
  std::string tag_name("Tracer");
  rval = mb.tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT);
  CHECK_ERR(rval);

  // tagElem is the average computed at each element, from nodal values
  Tag tagElem = 0;
  std::string tag_name2("TracerAverage");
  rval = mb.tag_get_handle(tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT);
  CHECK_ERR(rval);

  // area of the euler element is fixed, store it; it is used to recompute the averages at each
  // time step
  Tag tagArea = 0;
  std::string tag_name4("Area");
  rval = mb.tag_get_handle(tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT);
  CHECK_ERR(rval);

  // add a field value, quasi smooth first
  rval = add_field_value(&mb, euler_set, rank, tagTracer, tagElem, tagArea);
  CHECK_ERR(rval);

  // iniVals are used for 1-norm error computation
  Range redEls;
//.........这里部分代码省略.........

// Helper functions
void read_one_cell_var(bool rcbzoltan)
{
  Core moab;
  Interface& mb = moab;

  read_options = "PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL;NO_EDGES;VARIABLE=vorticity";
  if (rcbzoltan)
    read_options = "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;NO_EDGES;VARIABLE=vorticity;DEBUG_IO=1";

  ErrorCode rval = mb.load_file(example, NULL, read_options.c_str());
  CHECK_ERR(rval);

  // Get local edges
  Range local_edges;
  rval = mb.get_entities_by_type(0, MBEDGE, local_edges);
  CHECK_ERR(rval);
  CHECK_EQUAL((size_t)0, local_edges.size());

  // Get local cells
  Range local_cells;
  rval = mb.get_entities_by_type(0, MBPOLYGON, local_cells);
  CHECK_ERR(rval);
  // No mixed elements
  CHECK_EQUAL((size_t)1, local_cells.psize());

  Tag gid_tag;
  rval = mb.tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gid_tag, MB_TAG_DENSE);
  CHECK_ERR(rval);

  std::vector<int> gids(local_cells.size());
  rval = mb.tag_get_data(gid_tag, local_cells, &gids[0]);
  Range local_cell_gids;
  std::copy(gids.rbegin(), gids.rend(), range_inserter(local_cell_gids));

  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
  int rank = pcomm->proc_config().proc_rank();

  // Make check runs this test on two processors
  if (2 == procs) {
    // Check tag for cell variable vorticity at timestep 0
    Tag vorticity_tag0;
    rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0);
    CHECK_ERR(rval);

    // Check tag for cell variable vorticity at timestep 1
    Tag vorticity_tag1;
    rval = mb.tag_get_handle("vorticity1", layers, MB_TYPE_DOUBLE, vorticity_tag1);
    CHECK_ERR(rval);

    // Get vorticity0 and vorticity1 tag values on 3 local cells
    double vorticity0_val[3 * layers];
    double vorticity1_val[3 * layers];

    if (rcbzoltan) {
      CHECK_EQUAL((size_t)14, local_cell_gids.psize());

      if (0 == rank) {
        CHECK_EQUAL((size_t)319, local_cells.size());
        CHECK_EQUAL((size_t)319, local_cell_gids.size());

        CHECK_EQUAL(3, (int)local_cell_gids[0]);
        CHECK_EQUAL(162, (int)local_cell_gids[159]);
        CHECK_EQUAL(642, (int)local_cell_gids[318]);

        EntityHandle cell_ents[] = {local_cells[0], local_cells[159], local_cells[318]};
        rval = mb.tag_get_data(vorticity_tag0, cell_ents, 3, vorticity0_val);
        CHECK_ERR(rval);

        // Timestep 0
        // Layer 0
        CHECK_REAL_EQUAL(-0.725999, vorticity0_val[0 * layers], eps);
        CHECK_REAL_EQUAL(-1.814997, vorticity0_val[1 * layers], eps);
        CHECK_REAL_EQUAL(-0.554888, vorticity0_val[2 * layers], eps);
        // Layer 1
        CHECK_REAL_EQUAL(-0.725989, vorticity0_val[0 * layers + 1], eps);
        CHECK_REAL_EQUAL(-1.814972, vorticity0_val[1 * layers + 1], eps);
        CHECK_REAL_EQUAL(-0.554881, vorticity0_val[2 * layers + 1], eps);

        rval = mb.tag_get_data(vorticity_tag1, cell_ents, 3, vorticity1_val);
        CHECK_ERR(rval);

        // Timestep 1
        // Layer 0
        CHECK_REAL_EQUAL(-0.706871, vorticity1_val[0 * layers], eps);
        CHECK_REAL_EQUAL(-1.767178, vorticity1_val[1 * layers], eps);
        CHECK_REAL_EQUAL(-0.540269, vorticity1_val[2 * layers], eps);
        // Layer 1
        CHECK_REAL_EQUAL(-0.706861, vorticity1_val[0 * layers + 1], eps);
        CHECK_REAL_EQUAL(-1.767153, vorticity1_val[1 * layers + 1], eps);
        CHECK_REAL_EQUAL(-0.540262, vorticity1_val[2 * layers + 1], eps);
      }
      else if (1 == rank) {
        CHECK_EQUAL((size_t)323, local_cells.size());
        CHECK_EQUAL((size_t)323, local_cell_gids.size());

        CHECK_EQUAL(1, (int)local_cell_gids[0]);
        CHECK_EQUAL(365, (int)local_cell_gids[161]);
        CHECK_EQUAL(557, (int)local_cell_gids[322]);
//.........这里部分代码省略.........

void test_read_eul_onevar() 
{
  Core moab;
  Interface& mb = moab;
  std::string opts;
  ErrorCode rval = get_options(opts);
  CHECK_ERR(rval);

  opts += std::string(";VARIABLE=T");
  rval = mb.load_file(example_eul, NULL, opts.c_str());
  CHECK_ERR(rval);

  // Check for proper tags
  Tag Ttag0, Ttag1;
  rval = mb.tag_get_handle("T0", levels, MB_TYPE_DOUBLE, Ttag0);
  CHECK_ERR(rval);

  rval = mb.tag_get_handle("T1", levels, MB_TYPE_DOUBLE, Ttag1);
  CHECK_ERR(rval);

  // Check values of tag T0 (first level) at some strategically chosen places below
#ifdef MOAB_HAVE_MPI
  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int rank = pcomm->proc_config().proc_rank();
  int procs = pcomm->proc_config().proc_size();
#else
  int rank = 0;
  int procs = 1;
#endif

  const double eps = 0.0001;
  double val[8 * levels];

  if (1 == procs) {
    Range global_quads;
    rval = mb.get_entities_by_type(0, MBQUAD, global_quads);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)4608, global_quads.size());

    EntityHandle gloabl_quad_ents[] = {global_quads[0], global_quads[2255], global_quads[2304], global_quads[4559],
                                       global_quads[48], global_quads[2303], global_quads[2352], global_quads[4607]};
    rval = mb.tag_get_data(Ttag0, &gloabl_quad_ents[0], 8, val);

    CHECK_REAL_EQUAL(252.8529, val[0 * levels], eps); // First global quad
    CHECK_REAL_EQUAL(234.8390, val[1 * levels], eps); // 2256th global quad
    CHECK_REAL_EQUAL(232.6458, val[2 * levels], eps); // 2305th global quad
    CHECK_REAL_EQUAL(205.3905, val[3 * levels], eps); // 4560th global quad
    CHECK_REAL_EQUAL(252.7116, val[4 * levels], eps); // 49th global quad
    CHECK_REAL_EQUAL(232.6670, val[5 * levels], eps); // 2304th global quad
    CHECK_REAL_EQUAL(234.6922, val[6 * levels], eps); // 2353th global quad
    CHECK_REAL_EQUAL(200.6828, val[7 * levels], eps); // Last global quad
  }
  else if (2 == procs) {
    Range local_quads;
    rval = mb.get_entities_by_type(0, MBQUAD, local_quads);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)2304, local_quads.size());

    EntityHandle local_quad_ents[] = {local_quads[0], local_quads[1151], local_quads[1152], local_quads[2303]};
    rval = mb.tag_get_data(Ttag0, &local_quad_ents[0], 4, val);

    if (0 == rank) {
      CHECK_REAL_EQUAL(252.8529, val[0 * levels], eps); // First local quad, first global quad
      CHECK_REAL_EQUAL(234.8390, val[1 * levels], eps); // Median local quad, 2256th global quad
      CHECK_REAL_EQUAL(232.6458, val[2 * levels], eps); // Median local quad, 2305th global quad
      CHECK_REAL_EQUAL(205.3905, val[3 * levels], eps); // Last local quad, 4560th global quad
    }
    else if (1 == rank) {
      CHECK_REAL_EQUAL(252.7116, val[0 * levels], eps); // First local quad, 49th global quad
      CHECK_REAL_EQUAL(232.6670, val[1 * levels], eps); // Median local quad, 2304th global quad
      CHECK_REAL_EQUAL(234.6922, val[2 * levels], eps); // Median local quad, 2353th global quad
      CHECK_REAL_EQUAL(200.6828, val[3 * levels], eps); // Last local quad, last global quad
    }
  }
}

void test_read_all()
{
  Core moab;
  Interface& mb = moab;

  std::string opts;
  get_options(opts);

  // Read mesh and read all variables at all timesteps
  ErrorCode rval = mb.load_file(example, 0, opts.c_str());
  CHECK_ERR(rval);

#ifdef MOAB_HAVE_MPI
  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
#else
  int procs = 1;
#endif

  // Make check runs this test on one processor
  if (1 == procs) {
    // For u, wind and vorticity, check tag values on two entities
    double val[2 * layers];

    // Check tags for vertex variable u
    Tag u_tag0, u_tag1;
    rval = mb.tag_get_handle("u0", layers, MB_TYPE_DOUBLE, u_tag0);
    CHECK_ERR(rval);
    rval = mb.tag_get_handle("u1", layers, MB_TYPE_DOUBLE, u_tag1);
    CHECK_ERR(rval);

    // Get vertices (1280 vertices)
    Range verts;
    rval = mb.get_entities_by_type(0, MBVERTEX, verts);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)1280, verts.size());
    CHECK_EQUAL((size_t)1, verts.psize());

    // Check u tag values on first and last vertices
    EntityHandle vert_ents[] = {verts[0], verts[1279]};

    // Only check first two layers
    // Timestep 0
    rval = mb.tag_get_data(u_tag0, vert_ents, 2, val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(-4.839992, val[0 * layers], eps);
    CHECK_REAL_EQUAL(-3.699257, val[1 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(-4.839925, val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(-3.699206, val[1 * layers + 1], eps);

    // Timestep 1
    rval = mb.tag_get_data(u_tag1, vert_ents, 2, val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(-4.712473, val[0 * layers], eps);
    CHECK_REAL_EQUAL(-3.601793, val[1 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(-4.712409, val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(-3.601743, val[1 * layers + 1], eps);

    // Check tags for edge variable wind
    Tag wind_tag0, wind_tag1;
    rval = mb.tag_get_handle("wind0", layers, MB_TYPE_DOUBLE, wind_tag0);
    CHECK_ERR(rval);
    rval = mb.tag_get_handle("wind1", layers, MB_TYPE_DOUBLE, wind_tag1);
    CHECK_ERR(rval);

    // Get edges (1920 edges)
    Range edges;
    rval = mb.get_entities_by_type(0, MBEDGE, edges);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)1920, edges.size());
    CHECK_EQUAL((size_t)1, edges.psize());

    // Check wind tag values on first and last edges
    EntityHandle edge_ents[] = {edges[0], edges[1919]};

    // Only check first two layers
    // Timestep 0
    rval = mb.tag_get_data(wind_tag0, edge_ents, 2, val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(-5.081991, val[0 * layers], eps);
    CHECK_REAL_EQUAL(-6.420274, val[1 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(-5.081781, val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(-6.419831, val[1 * layers + 1], eps);

    // Timestep 1
    rval = mb.tag_get_data(wind_tag1, edge_ents, 2, val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(-4.948097, val[0 * layers], eps);
    CHECK_REAL_EQUAL(-6.251121, val[1 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(-4.947892, val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(-6.250690, val[1 * layers + 1], eps);

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

void test_read_onevar()
{
  Core moab;
  Interface& mb = moab;

  std::string opts;
  get_options(opts);

  // Read mesh and read cell variable vorticity at all timesteps
  opts += ";VARIABLE=vorticity";
  ErrorCode rval = mb.load_file(example, NULL, opts.c_str());
  CHECK_ERR(rval);

#ifdef MOAB_HAVE_MPI
  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
#else
  int procs = 1;
#endif

  // Make check runs this test on one processor
  if (1 == procs) {
    // Check vorticity tags
    Tag vorticity_tag0, vorticity_tag1;
    rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0);
    CHECK_ERR(rval);
    rval = mb.tag_get_handle("vorticity1", layers, MB_TYPE_DOUBLE, vorticity_tag1);
    CHECK_ERR(rval);

    // Get cells (12 pentagons and 630 hexagons)
    Range cells;
    rval = mb.get_entities_by_type(0, MBPOLYGON, cells);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)642, cells.size());

    // GCRM pentagons are always padded to hexagons
    CHECK_EQUAL((size_t)1, cells.psize());

    // Check vorticity tag values on 4 cells: first cell, two median cells, and last cell
    EntityHandle cell_ents[] = {cells[0], cells[320], cells[321], cells[641]};
    double vorticity_val[4 * layers];

    // Only check first two layers
    // Timestep 0
    rval = mb.tag_get_data(vorticity_tag0, cell_ents, 4, vorticity_val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(3.629994, vorticity_val[0 * layers], eps);
    CHECK_REAL_EQUAL(0.131688, vorticity_val[1 * layers], eps);
    CHECK_REAL_EQUAL(-0.554888, vorticity_val[2 * layers], eps);
    CHECK_REAL_EQUAL(-0.554888, vorticity_val[3 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(3.629944, vorticity_val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(0.131686, vorticity_val[1 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.554881, vorticity_val[2 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.554881, vorticity_val[3 * layers + 1], eps);

    // Timestep 1
    rval = mb.tag_get_data(vorticity_tag1, cell_ents, 4, vorticity_val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(3.534355, vorticity_val[0 * layers], eps);
    CHECK_REAL_EQUAL(0.128218, vorticity_val[1 * layers], eps);
    CHECK_REAL_EQUAL(-0.540269, vorticity_val[2 * layers], eps);
    CHECK_REAL_EQUAL(-0.540269, vorticity_val[3 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(3.534306, vorticity_val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(0.128216, vorticity_val[1 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.540262, vorticity_val[2 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.540262, vorticity_val[3 * layers + 1], eps);
  }
}

void test_read_novars()
{
  Core moab;
  Interface& mb = moab;

  // Need a file set for nomesh to work right
  EntityHandle file_set;
  ErrorCode rval = mb.create_meshset(MESHSET_SET, file_set);
  CHECK_ERR(rval);

  std::string orig, opts;
  get_options(orig);
  CHECK_ERR(rval);

  // Read header info only, no mesh, no variables
  opts = orig + ";NOMESH;VARIABLE=";
  rval = mb.load_file(example, &file_set, opts.c_str());
  CHECK_ERR(rval);

  // Read mesh, but still no variables
  opts = orig + ";VARIABLE=";
  rval = mb.load_file(example, &file_set, opts.c_str());
  CHECK_ERR(rval);

  // Check vorticity tags
  Tag vorticity_tag0, vorticity_tag1;
  rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0);
  // Tag vorticity0 should not exist
  CHECK_EQUAL(rval, MB_TAG_NOT_FOUND);
  rval = mb.tag_get_handle("vorticity1", layers, MB_TYPE_DOUBLE, vorticity_tag1);
  // Tag vorticity1 should not exist
  CHECK_EQUAL(rval, MB_TAG_NOT_FOUND);

  // Read vorticity at 1st timestep, no need to read mesh
  opts = orig + ";VARIABLE=vorticity;TIMESTEP=0;NOMESH";
  rval = mb.load_file(example, &file_set, opts.c_str());
  CHECK_ERR(rval);

  // Check vorticity tags again
  rval = mb.tag_get_handle("vorticity0", layers, MB_TYPE_DOUBLE, vorticity_tag0);
  // Tag vorticity0 should exist at this time
  CHECK_ERR(rval);
  // Tag vorticity1 should still not exist
  rval = mb.tag_get_handle("vorticity1", layers, MB_TYPE_DOUBLE, vorticity_tag1);
  CHECK_EQUAL(rval, MB_TAG_NOT_FOUND);

  // Read vorticity at 2nd timestep, no need to read mesh
  opts = orig + ";VARIABLE=vorticity;TIMESTEP=1;NOMESH";
  rval = mb.load_file(example, &file_set, opts.c_str());
  CHECK_ERR(rval);

  // Check tag vorticity1 again
  rval = mb.tag_get_handle("vorticity1", layers, MB_TYPE_DOUBLE, vorticity_tag1);
  // Tag vorticity1 should exist at this time
  CHECK_ERR(rval);

#ifdef MOAB_HAVE_MPI
  ParallelComm* pcomm = ParallelComm::get_pcomm(&mb, 0);
  int procs = pcomm->proc_config().proc_size();
#else
  int procs = 1;
#endif

  // Make check runs this test on one processor
  if (1 == procs) {
    // Get cells (12 pentagons and 630 hexagons)
    Range cells;
    rval = mb.get_entities_by_type(file_set, MBPOLYGON, cells);
    CHECK_ERR(rval);
    CHECK_EQUAL((size_t)642, cells.size());

    // GCRM pentagons are always padded to hexagons
    CHECK_EQUAL((size_t)1, cells.psize());

    // Check vorticity tag values on 4 cells: first cell, two median cells, and last cell
    EntityHandle cell_ents[] = {cells[0], cells[320], cells[321], cells[641]};
    double vorticity_val[4 * layers];

    // Only check first two layers
    // Timestep 0
    rval = mb.tag_get_data(vorticity_tag0, cell_ents, 4, vorticity_val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(3.629994, vorticity_val[0 * layers], eps);
    CHECK_REAL_EQUAL(0.131688, vorticity_val[1 * layers], eps);
    CHECK_REAL_EQUAL(-0.554888, vorticity_val[2 * layers], eps);
    CHECK_REAL_EQUAL(-0.554888, vorticity_val[3 * layers], eps);
    // Layer 1
    CHECK_REAL_EQUAL(3.629944, vorticity_val[0 * layers + 1], eps);
    CHECK_REAL_EQUAL(0.131686, vorticity_val[1 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.554881, vorticity_val[2 * layers + 1], eps);
    CHECK_REAL_EQUAL(-0.554881, vorticity_val[3 * layers + 1], eps);

    // Timestep 1
    rval = mb.tag_get_data(vorticity_tag1, cell_ents, 4, vorticity_val);
    CHECK_ERR(rval);
    // Layer 0
    CHECK_REAL_EQUAL(3.534355, vorticity_val[0 * layers], eps);
    CHECK_REAL_EQUAL(0.128218, vorticity_val[1 * layers], eps);
    CHECK_REAL_EQUAL(-0.540269, vorticity_val[2 * layers], eps);
//.........这里部分代码省略.........