sptlrgn.cpp

pic 模拟程序！面向对象

       //      strstream s;
#ifdef HAVE_SSTREAM
      stringstream s;
#else
      strstream s;
#endif
    
      // s << spGroupName <<"/pGrp"<< ptclGrpCount <<ends;
      s << speciesIter.current()->get_name() << "/pGrp"<< ptclGrpCount <<ends;
      string pGroupName = s.str();
           
      pgscan.current()->dumpH5(dumpObj, grid,pGroupName);
      ptclGrpCount++;
    }
    
    attrdataspaceId = H5Screate(H5S_SCALAR);
    attributeId = H5Acreate(spGroupId, "nPtclGrp",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
    status = H5Awrite(attributeId, H5T_NATIVE_INT, &ptclGrpCount );
    status = H5Aclose(attributeId);
    status = H5Sclose(attrdataspaceId);

    H5Gclose(spGroupId);
    isp++;
  }
//  H5Dclose(datasetId);
  H5Gclose(groupId);
  H5Fclose(fileId);
   
  return (1);
}
#endif

int SpatialRegion::Dump(FILE* DMPFile)
{

  cerr << "Entering sptlrgn::Dump"<<endl;

  static int end=-1;
  int nn;
  FILE *DropFile;
  // char dropfile[50]; // this is too small
  char dropfile[180];
  XGWrite(&frandseed, 4, 1, DMPFile, "long");
	
  int RnameLength = Name.length();
  // Write the length of the region name
  XGWrite(&RnameLength,4,1,DMPFile,"int");

  // Write the region name
  XGWrite(Name.c_str(),1,RnameLength+1,DMPFile,"char");

  // Write simulation time
  XGWrite(&t, ScalarInt, 1, DMPFile, ScalarChar);
	 
  // Write the physical dimensions of the system
  {
    Scalar xw[4];  // in this format:  xs, ys, xf, yf
    Vector2 **X=grid->getX();
    xw[0] = X[0][0].e1();
    xw[1] = X[0][0].e2();
    xw[2] = X[getJ()][getK()].e1();
    xw[3] = X[getJ()][getK()].e2();
    XGWrite(xw,ScalarInt,4,DMPFile,ScalarChar);
  }



  // make a drop file entry
  nn = strlen(theDumpFile);
  strcpy(dropfile, theDumpFile);

  dropfile[nn-2] = 'r';
  if ((DropFile = fopen(dropfile, "a")) == NULL)
    {
      printf("filename = %s\n", dropfile);
      puts("Drop: open failed error#1");
      return (0);
    }
  fprintf(DropFile, "%g ", t);

  for (int i=0; i<nSpecies; i++){
    int number = CountParticles(i);
    fprintf(DropFile, " %d ", number);
  }
  
  // Drop the number of shifts the spatial region has done.
  fprintf(DropFile," %lu ", shiftNum);

  fprintf(DropFile, "\n");
  fclose(DropFile);
	
 
  //  Dump Boundaries
  
  // dump the # of bounaries

  { int n;
  n = boundaryList->nItems();


  XGWrite(&n,4,1,DMPFile,"int");

  }
  oopicListIter<Boundary>	nextb(*boundaryList);
  int nnnn = 0;
  for (nextb.restart(); !nextb.Done(); nextb++)
    {
      if (!(nextb.current()->Dump(DMPFile)))
        printf("Dump: dump failed to save boundaries");
    }
	
  //  Dump asked for diagnostic
  // dump the # of diagnostics
  { int n;
  n = diagList->nItems();
  XGWrite(&n,4,1,DMPFile,"int");
  }
  oopicListIter<Diag>	nextd(*diagList);
  for (nextd.restart(); !nextd.Done(); nextd++)
    {
      if (!(nextd.current()->Dump(DMPFile)))
        cout << "Dump: dump failed to save diagnostics." << endl;
 }	

  end = -1;
  
  fields->Dump(DMPFile);


  //
  // dump of the NGDs
  // 
  oopicListIter<NGD> NGDIter(NGDList);

  for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ){ 
    NGDIter.current()->Dump(DMPFile); 
  }
  XGWrite(&nSpecies, 4, 1, DMPFile, "int");

  for (int isp=0; isp<nSpecies; isp++){
    oopicListIter<ParticleGroup> pgscan(*particleGroupList[isp]);
    for(pgscan.restart();!pgscan.Done();pgscan++)
      pgscan.current()->Dump(DMPFile,grid);
  }

  XGWrite(&end, 4, 1, DMPFile, "int");  //end of particle marker
  XGWrite(&end, 4, 1, DMPFile, "int");
  XGWrite(&end, 4, 1, DMPFile, "int");	
 
  return (1);
	
}


int SpatialRegion::Restore_2_00(FILE* DMPFile)
  throw(Oops){
  int j1test, k1test, j2test, k2test;
  int found;

  XGRead(&frandseed, 4, 1, DMPFile, "long");
	
  cerr << "\nentered  SpatialRegion::Restore_2_00\n";

  int ReadNameLength;
  char buf[512];
	
  XGRead(&ReadNameLength,4,1,DMPFile,"int");   // Read length of the region name
  XGRead(buf,1,ReadNameLength+1,DMPFile,"char"); // Read region name

  XGRead(&t, ScalarInt, 1, DMPFile, ScalarChar); // Read simulation time

  //  Restore Boundaries must be in the same order as dumped
  // Should give them titles as in diagnostics

  XGRead(&j1test, 4, 1, DMPFile, "int");
  XGRead(&k1test, 4, 1, DMPFile, "int");
  XGRead(&j2test, 4, 1, DMPFile, "int");
  XGRead(&k2test, 4, 1, DMPFile, "int");  
  oopicListIter<Boundary>	nextb(*boundaryList);
  for (nextb.restart(); !nextb.Done(); nextb++){
    if (!(nextb.current()->Restore_2_00(DMPFile, j1test, k1test, j2test, k2test))){
      stringstream ss (stringstream::in | stringstream::out);
      ss << "SpatialRegion::Restore_2_00: Error: \n"<<
        "Input file has boundaries that can't be matched to dump file." << endl;
     
      std::string msg;
      ss >> msg;
      Oops oops(msg);
      throw oops;    // exit() dump::Restore  dump::Restore_2_00

    }
    XGRead(&j1test, 4, 1, DMPFile, "int");
    XGRead(&k1test, 4, 1, DMPFile, "int");
    XGRead(&j2test, 4, 1, DMPFile, "int");
    XGRead(&k2test, 4, 1, DMPFile, "int");
  } 

	
  //  Restore diagnostic
  XGRead(&j1test, 4, 1, DMPFile, "int");
  XGRead(&k1test, 4, 1, DMPFile, "int");
  XGRead(&j2test, 4, 1, DMPFile, "int");
  XGRead(&k2test, 4, 1, DMPFile, "int"); 
  int titleLength;
  XGRead(&titleLength, 4,1,DMPFile,"int");
  XGRead(buf,1,titleLength+1, DMPFile, "char");
  ostring titleTest=ostring(buf);
  while (j1test>=0){
    found = FALSE;
    oopicListIter<Diag>	nextd(*diagList);
    for (nextd.restart(); !nextd.Done(); nextd++){
      if (nextd.current()->Restore_2_00(DMPFile, j1test, k1test, j2test, k2test, titleTest)){
        found = TRUE;
        break;
      }
    }
    if(!found){
       printf("Input file has diagnostics that can't be matched to dump file.");
      sprintf(buf, "Assuming that you dropped %s \n",titleTest());
      sprintf(buf, "If you 'save' again %s's data will be lost. \n", titleTest());
      //			cout << "Assuming that you dropped " << (titleTest.()) << endl;
      //			cout << "If you 'save' again " << (titleTest.()) << "'s data will be lost. << endl;
    }	
		
    XGRead(&j1test, 4, 1, DMPFile, "int");
    XGRead(&k1test, 4, 1, DMPFile, "int");
    XGRead(&j2test, 4, 1, DMPFile, "int");
    XGRead(&k2test, 4, 1, DMPFile, "int"); 
    XGRead(&titleLength, 4,1,DMPFile,"int");
    XGRead(buf,1,titleLength+1, DMPFile, "char");
    titleTest=ostring(buf);
  }
 
  fields->Restore_2_00(DMPFile);

  int speciesID;
  Scalar np2c0;
  int vary_np2c;
  Species* species=0;

  XGRead(&nSpecies, 4, 1, DMPFile, "int");
  XGRead(&speciesID, 4, 1, DMPFile, "int");
  while (!(speciesID==-1)){

    oopicListIter<Species> speciesIter(*speciesList);
    for (speciesIter.restart(); !speciesIter.Done(); speciesIter++)
      if (speciesIter.current()->getID() == speciesID)
        species = speciesIter.current();

    XGRead(&vary_np2c, 4, 1, DMPFile, "int");
    XGRead(&np2c0, ScalarInt, 1, DMPFile, ScalarChar);

    fprintf(stderr, "\nduplicateParticles=%d, np2c0=%f", 
	    duplicateParticles, np2c0);
		
    // Making a new particleGroup
		
    ParticleGroup* pg;
    

    if(getInfiniteBFlag() && getElectrostaticFlag())
      pg = new ParticleGroupIBES(maxN, species, np2c0, vary_np2c);
    else if (getInfiniteBFlag()) // relativistic Infinite B mover
      pg = new ParticleGroupIBEM(maxN, species, np2c0, vary_np2c);  
    else if (getElectrostaticFlag()) // add ES group
      pg = new ParticleGroupES(maxN, species, np2c0, vary_np2c);
    else if (getSynchRadiationFlag()) // add SR group
      pg = new ParticleGroupSR(maxN, species, np2c0, vary_np2c);
    else if(getNonRelFlag())
      pg = new ParticleGroupNR(maxN, species, np2c0, vary_np2c);
    else        //	add EM group
      pg = new ParticleGroup(maxN, species, np2c0, vary_np2c);

    pg->Restore_2_00(DMPFile);

    ParticleGroup** k;
    if(duplicateParticles > 0) 
      { 
	fprintf(stderr, "\nduplicateParticles = %d", duplicateParticles);
	k = pg->duplicateParticles(duplicateParticles); 
	if(k != NULL)
	  {
	    for(int i=0; i < duplicateParticles; i++)
	      {
		particleGroupList[speciesID]->add(k[i]);
	      }
	  }
	else
	  {
	    fprintf(stderr, "\nError duplicating particles");
	  }
      }
    particleGroupList[speciesID]->add(pg);
    
    XGRead(&speciesID, 4, 1, DMPFile, "int");
  }

  return (1);

}

#ifdef HAVE_HDF5
int SpatialRegion::restoreH5(dumpHDF5 restoreObj)
throw(Oops){

  // cerr << "entered SpatialRegion::restoreH5().\n"; 

  hid_t fileId, groupId; // datasetId,dataspaceId;
  herr_t status;
  string outFile;
//  hsize_t rank; 
//  hsize_t     dimsf[1];   
//  hsize_t dims;
  hid_t scalarType = restoreObj.getHDF5ScalarType();
  hid_t attrId;// attrspaceId
  hid_t bMetaGroupId; 
  string restoreFile = restoreObj.getDumpFileName();

  fileId = H5Fopen(restoreObj.getDumpFileName().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  if(fileId<0) {
    std::cerr << "Unable to open file '" << restoreObj.getDumpFileName()
	      << "'.  Error # " << fileId << std::endl;
    return(-1);
  }
  
// Open the group for the region
  groupId = H5Gopen(fileId, "/SpatialRegion");

// Read in the field data from the file

// Read the randomseed attribute
  attrId = H5Aopen_name(groupId, "frandseed");
  status = H5Aread(attrId,H5T_NATIVE_INT, &frandseed);
  //  cerr << "read random seed: "<<frandseed<<endl;
// Close the attribute
  status = H5Aclose(attrId);

// Read the time attribute
  attrId = H5Aopen_name(groupId, "simulationTime");
  status = H5Aread(attrId,scalarType , &t);
  // cerr << " read time: "<<time<<endl;
// Close the attribute
  status = H5Aclose(attrId);

  //  cerr << "read time.\n";

//read dimensions of region
  attrId = H5Aopen_name(groupId, "dims");
  // attrspaceId = H5Aget_space(attrId );
  Scalar xw[4];
  status = H5Aread(attrId, scalarType, xw);
  H5Aclose(attrId);
 
// Read the number boundaries attribute
  int numBoundaries;
  attrId = H5Aopen_name(groupId, "numBoundaries");
  status = H5Aread(attrId,H5T_NATIVE_INT , &numBoundaries);
  //  cerr << " read numBoundaries: "<<numBoundaries<<endl;
  status = H5Aclose(attrId);

  // status = H5Dclose(datasetId);
   
  bMetaGroupId = H5Gopen(fileId, "/Boundaries");
  for(int i=1;i<numBoundaries;i++) {
//    Scalar A1,A2,B1,B2;  // phys. dimensions of boundary
    Scalar bDims[4];
    hid_t bGroupId;
    int BoundaryType; // Type of this boundary
    int WasRestored=0; //Did one of the boundaries restore from this boundary?
    int nQuads;  // number of ordered quads of floats to restore
    
    // read in the boundary extent
   //read dimensions of boundary

    //    strstream s;
#ifdef HAVE_SSTREAM
    stringstream s;
#else
    strstream s;
#endif

    s <<  "boundary" << i <<ends;    // ends provides the null terminator (don`t forget that)
    string bname = s.str();

    bGroupId = H5Gopen(bMetaGroupId, bname.c_str());

    bname = "/Boundaries/" + bname;
    attrId = H5Aopen_name(bGroupId, "dims" );
    // attrspaceId = H5Dget_space(datasetId );
    status = H5Aread(attrId, scalarType, bDims);
    H5Aclose(attrId);
    
// read in boundary type
    attrId = H5Aopen_name(bGroupId, "boundaryType" );
    // attrspaceId = H5Dget_space(datasetId );
    status = H5Aread(attrId, H5T_NATIVE_INT, &BoundaryType);
    H5Aclose(attrId);
    
    attrId = H5Aopen_name(bGroupId, "nQuads" );
    // attrspaceId = H5Dget_space(datasetId );
    status = H5Aread(attrId, H5T_NATIVE_INT, &nQuads);
    H5Aclose(attrId);

    if(nQuads>0) {
// try to find a matching boundary
      oopicListIter<Boundary>	nextb(*boundaryList);
      for (nextb.restart(); !nextb.Done()&& WasRestored==0; nextb++){
	//if(BoundaryType != -4)
		WasRestored=nextb.current()->restoreH5(restoreObj,BoundaryType,bname,nQuads);
      
 }
    }
  H5Gclose(bGroupId);
  }
  H5Gclose(bMetaGroupId);

  //
// restore diagnostics (empty for now)
  //

  //
// restore fields
  //
   fields->restoreH5(restoreObj);

   //
// restore NGD
   //

  oopicListIter<NGD> NGDIter(NGDList);
  int nnnn = 0;
  for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ){ 
    //    strstream s;
#ifdef HAVE_SSTREAM
    stringstream s;
#else
    strstream s;
#endif

    s << "NGD" << nnnn <<ends;   
    string ngdGrpName = s.str();
    try{
      NGDIter.current()->restoreH5(restoreObj,ngdGrpName); 
    }
    catch(Oops& oops){
      oops.prepend(" SpatialRegion::restoreH5: Error: \n");//dumpHDF5::restoreH5
      throw oops;
    }
    nnnn++;
}

   //
// restore species
   //
  //  cerr << "restoring Species";


  int speciesID;
  Scalar np2c0;
  int vary_np2c;
  

// Read the nSpecies attribute
  attrId = H5Aopen_name(groupId, "nSpecies");
  status = H5Aread(attrId,H5T_NATIVE_INT , &nSpecies);
  status = H5Aclose(attrId);

  // for(int i=0;i<nSpecies;i++)


  oopicListIter<Species> speciesIter(*speciesList);
  for (speciesIter.restart(); !speciesIter.Done(); speciesIter++)
    {
      
//        strstream sp;
//        sp << "/species" << i <<ends;
      ostring spGrpName = speciesIter.current()->get_name();
      Species* species = speciesIter.current();

      hid_t spGroupId = H5Gopen(fileId,speciesIter.current()->get_name().c_str());
// Read the nPtclGrp attribute
       int nPtclGrp;
       attrId = H5Aopen_name(spGroupId, "nPtclGrp");
       status = H5Aread(attrId,H5T_NATIVE_INT , &nPtclGrp);