exitport.cpp

pic 模拟程序！面向对象

//	ExitPort::emit() simply deletes Particles in its stack.  May add some
//	diagnostics for particles collected in the future.
#if !defined __linux__ && !defined _WIN32
#pragma argsused
#endif
/*
ParticleList& ExitPort::emit(Scalar t,Scalar dt)
{
	while(!particleList.isEmpty())
		{
			Particle* p = particleList.pop();
			delete p;
		}
	return particleList;
}
*/
#ifdef HAVE_HDF5
int ExitPort::dumpH5(dumpHDF5 &dumpObj,int number)
{

  //cerr << "entering ExitPort::dumpH5() . \n";
//  Scalar temp;
//  Scalar x,y;
  int length;
  
  hid_t fileId, groupId, datasetId,dataspaceId,HspaceId,HsetId  ;
  herr_t status;
  hsize_t rank; 
//  hsize_t     dimsf[1];   
  hsize_t dims;
  hid_t scalarType = dumpObj.getHDF5ScalarType();
//  int         attrData;
  hid_t  attrdataspaceId,attributeId;
#ifdef NEW_H5S_SELECT_HYPERSLAB_IFC
  hsize_t start[2]; /* Start of hyperslab */
#else
  hssize_t start[2]; /* Start of hyperslab */
#endif
  hsize_t count[2];  /* Block count */

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

  s << "boundary" << number <<ends; 
  string portname = s.str();
   portname = "/Boundaries/" + portname;
// Open the hdf5 file with write access
  fileId = H5Fopen(dumpObj.getDumpFileName().c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
  groupId = H5Gcreate(fileId,portname.c_str(),0);
// get the data   
 

  Grid *grid = fields->get_grid();
// Write the physical dimensions of the boundary
  
  Scalar xw[4];  // in this format:  xs, ys, xf, yf
  Vector2 **X=grid->getX();
  xw[0] = X[j1][k1].e1();
  xw[1] = X[j1][k1].e2();
  xw[2] = X[j2][k2].e1();
  xw[3] = X[j2][k2].e2();
  
/* Create the data space for the attribute. */
  dims = 4;
  attrdataspaceId = H5Screate_simple(1, &dims, NULL);
/* Create a dataset attribute. */
  attributeId = H5Acreate(groupId, "dims",scalarType, attrdataspaceId, H5P_DEFAULT);
/* Write the attribute data. */
  status = H5Awrite(attributeId, scalarType, xw);
/* Close the attribute. */
  status = H5Aclose(attributeId);
/* Close the dataspace. */
  status = H5Sclose(attrdataspaceId);

// write Boundary Type as attribute
  dims = 1;
  attrdataspaceId = H5Screate_simple(1, &dims, NULL);
  attributeId = H5Acreate(groupId, "boundaryType", H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, H5T_NATIVE_INT, &BoundaryType);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);

// write length as attribute --------
  length = MAX( (abs(k2-k1)),(abs(j2-j1))) +1;
  dims = 1;
  attrdataspaceId = H5Screate_simple(1, &dims, NULL);
  attributeId = H5Acreate(groupId, "nQuads", H5T_NATIVE_INT , attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, H5T_NATIVE_INT, &length);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);


  length = MAX( (abs(k2-k1)),(abs(j2-j1))) +1;
  rank = 2;
  hsize_t     dimsf2[2];
  dimsf2[0] = length;
  dimsf2[1] = 2;

  dataspaceId = H5Screate_simple(rank, dimsf2, NULL); 
  string datasetName = "Points";  
  datasetId = H5Dcreate(groupId, datasetName.c_str(), scalarType, dataspaceId,
			H5P_DEFAULT);
  
  HspaceId = H5Screate_simple(rank, dimsf2, NULL); 
  HsetId = H5Dcreate(groupId, "oldH", scalarType, HspaceId,
			H5P_DEFAULT);
  hsize_t dim2 = 2;
  hid_t slabId = H5Screate_simple(1, &dim2, NULL);

  if (alongx2())	// vertical
    {
      int kmax = abs(k2-k1) + 1;
      for (int k = 0; k<kmax; k++)
	{
	  Scalar data[2] = {X[j1][k].e1(),X[j1][k].e2()};  
	  start[0]  = 0;  count[0]  = 2;
	  status = H5Sselect_hyperslab(slabId, H5S_SELECT_SET, start, NULL, count, NULL);
	  start[0]  = k; start[1]  = 0;
	  count[0]  = 1; count[1]  = 2;     
	  status = H5Sselect_hyperslab(dataspaceId, H5S_SELECT_SET, start, NULL, count, NULL);
	  status = H5Dwrite(datasetId, scalarType, slabId, dataspaceId, H5P_DEFAULT, data);
     
	  data[0] = oldHonBoundary[k].e1();  
	  data[1] = oldHonBoundary[k].e2();   
	  start[0]  = 0;  count[0]  = 2;
	  status = H5Sselect_hyperslab(slabId, H5S_SELECT_SET, start, NULL, count, NULL);
	  start[0]  = k; start[1]  = 0;
	  count[0]  = 1; count[1]  = 2;     
	  status = H5Sselect_hyperslab(HspaceId, H5S_SELECT_SET, start, NULL, count, NULL);
	  status = H5Dwrite(HsetId, scalarType, slabId, HspaceId, H5P_DEFAULT, data);
	}
    }
  else   // horizontal
    {
      int jmax = abs(j2-j1) + 1;
      for (int j = 0; j<jmax; j++)
	{
	  Scalar data[2] = {X[j][k1].e1(),X[j][k1].e2()};  
	  start[0]  = 0;  count[0]  = 2;
	  status = H5Sselect_hyperslab(slabId, H5S_SELECT_SET, start, NULL, count, NULL);
	  start[0]  = j; start[1]  = 0;
	  count[0]  = 1; count[1]  = 2;     
	  status = H5Sselect_hyperslab(dataspaceId, H5S_SELECT_SET, start, NULL, count, NULL);
	  status = H5Dwrite(datasetId, scalarType, slabId, dataspaceId, H5P_DEFAULT, data);

	  data[0] = oldHonBoundary[j].e1();  
	  data[1] = oldHonBoundary[j].e2();   
	  start[0]  = 0;  count[0]  = 2;
	  status = H5Sselect_hyperslab(slabId, H5S_SELECT_SET, start, NULL, count, NULL);
	  start[0]  = j; start[1]  = 0;
	  count[0]  = 1; count[1]  = 2;     
	  status = H5Sselect_hyperslab(HspaceId, H5S_SELECT_SET, start, NULL, count, NULL);
	  status = H5Dwrite(HsetId, scalarType, slabId, HspaceId, H5P_DEFAULT, data);
	}
    }	
   
  H5Sclose(dataspaceId);
  H5Dclose(datasetId);   
  H5Sclose(HspaceId);
  H5Dclose(HsetId);   

  // write fields
   
  H5Gclose(groupId);
  H5Fclose(fileId);
  return(1);
}
#endif

int ExitPort::Dump(FILE * DMPFile)
{
  Scalar temp;
  Scalar x,y;
  int length;

	
  //#ifdef UNIX
  {
    Scalar xw[4];  // in this format:  xs, ys, xf, yf
    Vector2 **X=grid->getX();
    xw[0] = X[j1][k1].e1();
    xw[1] = X[j1][k1].e2();
    xw[2] = X[j2][k2].e1();
    xw[3] = X[j2][k2].e2();
    XGWrite(xw,ScalarInt,4,DMPFile,ScalarChar);
  }
  // write the BoundaryType
  XGWrite(&BoundaryType, 4, 1, DMPFile, "int");


  Vector2** X=fields->get_grid()->getX();
	
  length = MAX( (abs(k2-k1)),(abs(j2-j1))) +1;
  XGWrite(&length,4,1,DMPFile,"int");
  XGWrite(&length,4,1,DMPFile,"int");  // redundant, but necessary
  XGWrite(&length,4,1,DMPFile,"int");  // redundant, but necessary
  if (alongx2())	// vertical
    {
      int kmax = abs(k2-k1) + 1;
      for (int k = 0; k<kmax; k++)
        {
          x = X[j1][k].e1();
          y = X[j1][k].e2();
          XGWrite(&x,ScalarInt,1,DMPFile,ScalarChar);
          XGWrite(&y,ScalarInt,1,DMPFile,ScalarChar);
          temp = oldHonBoundary[k].e1();
          XGWrite(&temp, ScalarInt, 1, DMPFile, ScalarChar);
          temp = oldHonBoundary[k].e2();
          XGWrite(&temp, ScalarInt, 1, DMPFile, ScalarChar);
        }
    }
  else   // horizontal
    {
      int jmax = abs(j2-j1) + 1;
      for (int j = 0; j<jmax; j++)
        {
          x = X[j][k1].e1();
          y = X[j][k1].e2();
          XGWrite(&x,ScalarInt,1,DMPFile,ScalarChar);
          XGWrite(&y,ScalarInt,1,DMPFile,ScalarChar);
          temp = oldHonBoundary[j].e1();
          XGWrite(&temp, ScalarInt, 1, DMPFile, ScalarChar);
          temp = oldHonBoundary[j].e2();
          XGWrite(&temp, ScalarInt, 1, DMPFile, ScalarChar);
        }
    }	
  //#endif
  return(TRUE);
}



#ifdef HAVE_HDF5
int ExitPort::restoreH5(dumpHDF5 &restoreObj,int bType,string whichBoundary,int nQuads)
{
  //cerr << "entering ExitPort::restoreH5().\n";
//  int dummy;
  int i;
  int jl,kl;
  Scalar bDims[4];
  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  attrspaceId,
  hid_t attrId;
#ifdef NEW_H5S_SELECT_HYPERSLAB_IFC
  hsize_t offset[2]; /* Offset of hyperslab */
#else
  hssize_t offset[2]; /* Offset of hyperslab */
#endif
  hsize_t count[2];  /* Block count */


  fileId = H5Fopen(restoreObj.getDumpFileName().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  groupId = H5Gopen(fileId, whichBoundary.c_str());

  //read dimensions of exitport
  attrId = H5Aopen_name(groupId, "dims" );
  status = H5Aread(attrId, scalarType, bDims);
  H5Aclose(attrId);

 
  

  if(BoundaryType != bType) return FALSE;
  if(!onMe(bDims[0],bDims[1],bDims[2],bDims[3])) return FALSE;
    
  // OK, this dielectric must be us.

  // Open the dataset for reading.  Determine shape and set field size.
  datasetId = H5Dopen(groupId,"Points" );
  dataspaceId = H5Dget_space(datasetId );
  int dataRank = H5Sget_simple_extent_ndims(dataspaceId);
  hsize_t *sizes = new hsize_t(dataRank);
  int res = H5Sget_simple_extent_dims(dataspaceId, sizes, NULL);
  hsize_t dim2 = sizes[1]; //should be 2
  hid_t posSlabId =  H5Screate_simple(1, &dim2, NULL);

  hid_t HsetId = H5Dopen(groupId,"oldH" );
  hid_t  HspaceId = H5Dget_space(HsetId );

  for(i=0;i<nQuads;i++) {
    Scalar x,y,H1,H2;
    Vector2 GridLoc;
 
    Scalar pos[2];   
    offset[0]  = 0;  count[0]  = 2;
    status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, offset, NULL, count, NULL);
    offset[0]  = i; offset[1]  = 0;;
    count[0]  = 1; count[1]  = 2;    
    status = H5Sselect_hyperslab(dataspaceId, H5S_SELECT_SET, offset, NULL, count, NULL);

    status = H5Dread(datasetId, scalarType, posSlabId, dataspaceId, H5P_DEFAULT, pos);
	      
    x = pos[0];
    y = pos[1];

    Scalar H[2];   
    //    status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, offset, NULL, count, NULL);
    status = H5Sselect_hyperslab(HspaceId, H5S_SELECT_SET, offset, NULL, count, NULL);
    status = H5Dread(HsetId, scalarType, posSlabId, HspaceId, H5P_DEFAULT, H);

    H1 = H[0];
    H2 = H[1];
    // Find the j/k value of the nearest
    // point on this boundary
    GridLoc = fields->get_grid()->getNearestGridCoords(Vector2(x,y));
    jl = (int) (GridLoc.e1() + 0.5);
    kl = (int) (GridLoc.e2() + 0.5);

    // find the point on this boundary this GridLoc corresponds to
    if (alongx2())	// vertical
      {
        int kmax = abs(k2-k1) +1;
        int ki = kl - MIN(k1,k2);
        if(ki <= kmax && kl >= 0)
          oldHonBoundary[ki]= Vector2(H1,H2);

      }
    else   // horizontal
      {
        int jmax = abs(j2-j1) + 1;
        int ji = jl = MIN(j1,j2);
        if(ji <= jmax && ji >=0)
          oldHonBoundary[ji]= Vector2(H1,H2);
      }	
    // #endif
  }


  H5Sclose(dataspaceId);
  H5Dclose(datasetId);
  H5Sclose(HspaceId);
  H5Dclose(HsetId);
  H5Sclose(posSlabId);
  H5Gclose(groupId);
  H5Fclose(fileId);

  delete sizes; 
  return(1);
}
#endif



int ExitPort::Restore(FILE *DMPFile,int _BoundaryType, Scalar A1, 
                      Scalar A2, Scalar B1, Scalar B2, int nQuads)
{
	
  int dummy;
  int i;
  int jl,kl;
  // #ifdef UNIX
  // Two checks to see if this data is loadable by this boundary.
  if(BoundaryType != _BoundaryType) return FALSE;
  if(!onMe(A1,A2,B1,B2)) return FALSE;
	

  // next two don't matter
  XGRead(&dummy,4,1,DMPFile,"int");
  XGRead(&dummy,4,1,DMPFile,"int");
	
//  Scalar temp1, temp2;

  for(i=0;i<nQuads;i++) {
    Scalar x,y,H1,H2;
    Vector2 GridLoc;
    XGRead(&x,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&y,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&H1,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&H2,ScalarInt,1,DMPFile,ScalarChar);
    // Find the j/k value of the nearest
    // point on this boundary
    GridLoc = fields->get_grid()->getNearestGridCoords(Vector2(x,y));
    jl = (int) (GridLoc.e1() + 0.5);
    kl = (int) (GridLoc.e2() + 0.5);

    // find the point on this boundary this GridLoc corresponds to
    if (alongx2())	// vertical
      {
        int kmax = abs(k2-k1) +1;
        int ki = kl - MIN(k1,k2);
        if(ki <= kmax && kl >= 0)
          oldHonBoundary[ki]= Vector2(H1,H2);

      }
    else   // horizontal
      {
        int jmax = abs(j2-j1) + 1;
        int ji = jl = MIN(j1,j2);
        if(ji <= jmax && ji >=0)
          oldHonBoundary[ji]= Vector2(H1,H2);
      }	
    // #endif
  }
  return(TRUE);

}

int ExitPort::Restore_2_00(FILE * DMPFile, int j1test, 
                           int k1test, int j2test, int k2test)
{
  if ((j1 == j1test)&&(k1 == k1test)&&(j2 == j2test)&&(k2 == k2test))
    {
      //#ifdef UNIX
      Scalar temp1, temp2;

      if (alongx2())	// vertical
        {
          int kmax = abs(k2-k1) + 1;
          for (int k = 0; k<kmax; k++)
            {
              XGRead(&temp1, ScalarInt, 1, DMPFile, ScalarChar);
              XGRead(&temp2, ScalarInt, 1, DMPFile, ScalarChar);
              oldHonBoundary[k] = Vector2(temp1,temp2);
            }
        }
      else   // horizontal
        {
          int jmax = abs(j2-j1) + 1;
          for (int j = 0; j<jmax; j++)
            {
              XGRead(&temp1, ScalarInt, 1, DMPFile, ScalarChar);
              XGRead(&temp2, ScalarInt, 1, DMPFile, ScalarChar);
              oldHonBoundary[j] = Vector2(temp1,temp2);
            }
        }	
      //#endif
      return(TRUE);
    }
  else
    return(FALSE);
}