dielectr.cpp

pic 模拟程序！面向对象

/*
   ====================================================================

   dielectr.cpp

   0.90	(PeterM 1-9-95) Creation.
   0.91  (JV/PM Mon Jul  3 ) Used bilinear weighting of Q.
   0.992 (JohnV 08-09-95) Add DielectricRegion
   0.993 (JohnV 09-02-95) Move secondary params into setSecondary().
   0.994 (JohnV 10-04-95) Add secondary emission for DielectricRegion, add
   capability for secondary coeff. > 1.
   1.001	(JohnV 05-07-96) Add toNodes() to handle improved interpolation
   at the boundary.
   1.1   (JohnV 10-09-96) Add reflection = fraction of particles reflected.
   1.11  (JohnV 12-16-96) Add refMaxE = Energy ceiling for reflection.
   1.12  (JohnV 02-07-96) Repaired potential problem for edge detection in 
   DielectricRegion::collect().  Also modified offset for secondaries.
   1.13  (JohnV 06-03-97) Added charge collection for secondaries.
   2.0   (JohnV 08-19-97) Added multiple secondary packages.
   2.1   (JohnV 06-15-98) Added Secondary::threshold.
   2.11  (JohnV 06-25-98) Repaired bug for reflected particles at corners.
   2.2   (JohnV 09-01-98) Moved Secondary to secondary.cpp.
   2.3   (JohnV 04-Mar-2002) Added transparency.

   ====================================================================
   */

#include "dielectr.h"
#include "fields.h"
#include "ptclgrp.h"
#include "misc.h"
#include "psolve.h"
#include "secondary.h"

#ifdef _MSC_VER
using std::cout;
using std::cerr;
using std::endl;
#endif
Dielectric::Dielectric(oopicList <LineSegment> *segments, Scalar _er, 
  int _QuseFlag, Scalar reflection, Scalar refMaxE) 
  : Boundary(segments)
{
  setQuseFlag(_QuseFlag);
  if (QuseFlag){
    nPoints = MAX(abs(j1-j2), abs(k1-k2)) + 1;
    Q = new Scalar[nPoints];
    memset(Q, 0, nPoints*sizeof(Scalar));
  }
  else {
    Q = NULL;
    nPoints = 0;
  }
  epsilon = _er/iEPS0;
  BoundaryType = DIELECTRIC;
  BCType = DIELECTRIC_BOUNDARY;
  if (alongx2()) {
    _delta.set_e1((j1==0)? 1e-6 : normal*1e-6*j1);
    _unitNormal.set_e1(normal);
  }
  else {
    _delta.set_e2((k1==0) ? 1e-6 : normal*1e-6*k1);
    _unitNormal.set_e2(normal);
  }
}

Dielectric::~Dielectric()
{
  if (Q) {
    delete [] Q;
    Q = 0;
  }
  secondaryList.deleteAll();
}

void Dielectric::setPassives()
{
}

//  This function sets the voltage for the use of the Poisson solve

void Dielectric::putCharge_and_Current(Scalar t)
{
   int j=0;
   int k=0;
   int ko=0;
   int jo=0;
   if (QuseFlag)
      {
         Scalar **rho=fields->getrho();
         Scalar **cellvolume = fields->get_grid()->get_halfCellVolumes();
         
         if (alongx2())								//	vertical
            for (k=ko=MIN(k1,k2); k<MAX(k1,k2); k++)
               rho[j1][k]+=Q[k-ko]/cellvolume[j1][k];
         else											//	horizontal
            for (jo=j=MIN(j1,j2); j<MAX(j1,j2); j++)
               rho[j][k1]+=Q[j-jo]/cellvolume[j][k1];
      }
}

void Dielectric::putCharge(void)
{
   int j,k,ko,jo;
   if (QuseFlag){
      Scalar **Charge=fields->getQ();
      
      if (alongx2())								//	vertical
         for (k=ko=MIN(k1,k2); k<MAX(k1,k2); k++)
            Charge[j1][k]+=Q[k-ko];
      else											//	horizontal
         for (jo=j=MIN(j1,j2); j<MAX(j1,j2); j++)
            Charge[j][k1]+=Q[j-jo];
   }
}

void Dielectric::setFields(Fields &f)
{
   Boundary::setFields(f);
}

void Dielectric::addSecondary(Secondary* newSecondary)
{
   if (newSecondary->checkSpecies()) secondaryList.add(newSecondary);
   newSecondary->setBoundaryPtr(this);
}

//--------------------------------------------------------------------
// Dielectric::collect() handles weighting to charge to boundary, and
//	also secondary emission.  Some fraction of particles (reflection)
// can be reflected back, when energy <= refMaxE.

void	Dielectric::collect(Particle& p, Vector3& dxMKS)
{
   int j=0;
   Scalar w=0.;
   if ((Scalar)frand() <= reflection && p.energy_eV_MKS() <= refMaxE){  // reflect particles
      Vector3 u = p.get_u();
      Vector2 x = p.get_x();
      if (alongx2()){
         u.set_e1(-u.e1()); // reverse u comp. normal to bc
         dxMKS.set_e1(-dxMKS.e1()); // reverse the movement
         if (u.e1() < 0) x.set_e1(x.e1()*.99999);
         else x.set_e1(x.e1()*1.00001);
         if (x.e2()==(int)x.e2()) x.set_e2(x.e2()*1.00001); // move off corners
         else if (x.e2()==(int)x.e2()+1) x.set_e2(x.e2()*0.99999);
      }
      else {
         u.set_e2(-u.e2());
         dxMKS.set_e2(-dxMKS.e2()); // reverse the movement
         if (u.e2() < 0) x.set_e2(x.e2()*.99999);
         else x.set_e2(x.e2()*1.00001);
         if (x.e1()==(int)x.e1()) x.set_e1(x.e1()*1.00001); // move off corners
         else if (x.e1()==(int)x.e1()+1) x.set_e1(x.e1()*0.99999);
      }
      // note that the dt used below should probably be the remaining fraction for the
      // move rather than the full dt.
      Boundary *bPtr = fields->translateAccumulate(x, dxMKS, p.get_q()/fields->get_dt());
      p.set_u(u);
      p.set_x(x);
      if (bPtr) bPtr->collect(p, dxMKS);
      else particleList.add(&p);
      return;
   }
   else if ((Scalar)frand() < transparency) {
     // note that the dt used below should probably be the remaining fraction for the
     // move rather than the full dt.
     Vector2 x = p.get_x();
     Boundary *bPtr = fields->translateAccumulate(x, dxMKS, p.get_q()/fields->get_dt());
     p.set_x(x);
     if (bPtr) bPtr->collect(p, dxMKS);
     else particleList.add(&p);
     return;
   } 
   if (QuseFlag) {
      if (alongx2()){	// vertical
         j = (int) p.get_x().e2();
         w = p.get_x().e2() - j;
         j -= k1;
      }
      else {
         j = (int) p.get_x().e1();
         w = p.get_x().e1() - j;
         j -= j1;
      }
      Q[j] += p.get_q()*(1-w);
      Q[MIN(j+1, nPoints-1)] += p.get_q()*w;
   }
   Boundary::collect(p, dxMKS); // do statistics;
   if (secondaryList.nItems()){
      p.set_x(p.get_x() + _delta); // perturb off edge;
      oopicListIter<Secondary> sIter(secondaryList);
      for (; !sIter.Done(); sIter++){
         int n = sIter()->generateSecondaries(p, particleList);
         if (QuseFlag){
            Scalar q = n*p.get_np2c()*sIter()->get_q();
            Q[j] -= q*(1-w);
            Q[MIN(j+1, nPoints-1)] -= q*w;
         }
      }
   }
   delete &p;
}
#ifdef HAVE_HDF5
int Dielectric::dumpH5(dumpHDF5 &dumpObj,int number)
{
  
  //cerr << "Entered Dielectric::dumpH5.\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 = dumpObj.getHDF5ScalarType();
  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;    // ends provides the null terminator (don`t forget that)
  string dname = s.str();
    
//  // Open the hdf5 file with write accessx
  fileId = H5Fopen(dumpObj.getDumpFileName().c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
  dname = "/Boundaries/" + dname;
  groupId = H5Gcreate(fileId,dname.c_str() ,0); 
  
  string datasetName = "dims";  
  
  rank = 1;
  dims = 4;
  attrdataspaceId = H5Screate_simple(rank, &dims, NULL); 
  attributeId = H5Acreate(groupId, datasetName.c_str(), scalarType, attrdataspaceId,
			  H5P_DEFAULT);
// 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();
  
  status = H5Awrite(attributeId, scalarType, xw);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);

// dump 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);

// dump n as attribute
  dims = 1;
  int n = 0;
  attrdataspaceId = H5Screate_simple(1, &dims, NULL);
  attributeId = H5Acreate(groupId, "nQuads",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, H5T_NATIVE_INT,&nPoints );
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);
  
//dump positions, charges 


  if(nPoints > 0){

  rank = 2;
  hsize_t     dimsf2[2];
  dimsf2[0] = nPoints;
  dimsf2[1] = 2;

  //  Scalar data[nPoints][2];

  dataspaceId = H5Screate_simple(rank, dimsf2, NULL); 
  datasetName = "positions";  
  datasetId = H5Dcreate(groupId, datasetName.c_str(), scalarType, dataspaceId,
			H5P_DEFAULT);
  hsize_t dim2 = 2;
  hid_t posSlabId = H5Screate_simple(1, &dim2, NULL);


  Scalar x,y;
  for(int i=0; i<nPoints; i++) {
    int ko = MIN(k1,k2);
    int jo = MIN(j1,j2);
    if(alongx2())   // Vertical
      {
	x = X[j1][i+ko].e1();
	y = X[j1][i+ko].e2();
	
      }
    else  // Horizontal
      {
	x = X[i+jo][k1].e1();
	y = X[i+jo][1].e2();
      }
    Scalar pos[2] = {x,y};  
    start[0]  = 0;  count[0]  = 2;
    status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, start, NULL, count, NULL);
    start[0]  = i; start[1]  = 0;
    count[0]  = 1; count[1]  = 2;    
    status = H5Sselect_hyperslab(dataspaceId, H5S_SELECT_SET, start, NULL, count, NULL);
    status = H5Dwrite(datasetId, scalarType, posSlabId, dataspaceId, H5P_DEFAULT, pos);
  }
  
  H5Sclose(posSlabId);
  H5Sclose(dataspaceId);
  H5Dclose(datasetId);
  
  // dump charges
  rank = 1;
  dimsf[0] = nPoints;
  datasetName = "Q";  
  dataspaceId = H5Screate_simple(rank, dimsf, NULL); 
  datasetId = H5Dcreate(groupId, datasetName.c_str(), scalarType, dataspaceId,
			H5P_DEFAULT);
  
  status = H5Dwrite(datasetId, scalarType, H5S_ALL, H5S_ALL,
		    H5P_DEFAULT, Q);
  H5Sclose(dataspaceId);
  H5Dclose(datasetId);
  }


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

int Dielectric::Dump(FILE * DMPFile)
{
  //#ifdef UNIX
  {
    Scalar xw[4];  // in this format:  xs, ys, xf, yf
    Vector2 **X=fields->get_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");
  
   XGWrite(&nPoints,4,1,DMPFile,"int");
   XGWrite(&nPoints,4,1,DMPFile,"int");  // redundant, but necessary
   XGWrite(&nPoints,4,1,DMPFile,"int");  // redundant, but necessary
   
   Vector2 **X = fields->get_grid()->getX();
   Scalar x,y;
   for(int i=0; i<nPoints; i++) {
     int ko = MIN(k1,k2);
     int jo = MIN(j1,j2);
     if(alongx2())   // Vertical
       {
	 x = X[j1][i+ko].e1();
	 y = X[j1][i+ko].e2();
       }
     else  // Horizontal
         {
	   x = X[i+jo][k1].e1();
	   y = X[i+jo][k1].e2();
         }
     XGWrite(&x,ScalarInt,1,DMPFile,ScalarChar);
     XGWrite(&y,ScalarInt,1,DMPFile,ScalarChar);
     XGWrite(&Q[i], ScalarInt, 1, DMPFile, ScalarChar);
     XGWrite(&Q[i], ScalarInt, 1, DMPFile, ScalarChar);  //redundant, but needed for restore
   }
   //#endif
   return(TRUE);
}


#ifdef HAVE_HDF5
int Dielectric::restoreH5(dumpHDF5 &restoreObj,int bType,string whichBoundary,int nQuads)
{
  //cerr << "entering Dielectric::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 attrId;// attrspaceId,
#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());

 
  
  attrId = H5Aopen_name(groupId, "dims" );
  status = H5Aread(attrId, scalarType, bDims);
  H5Aclose(attrId);
  
  if(BoundaryType != bType){
    H5Gclose(groupId);
    H5Fclose(fileId);
    return FALSE;
  }
  if(!onMe(bDims[0],bDims[1],bDims[2],bDims[3])) {
    H5Gclose(groupId);
    H5Fclose(fileId);
    return FALSE;
  }
  



//    // OK, this dielectric must be us.
  

  datasetId = H5Dopen(groupId,"positions");
  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);
  int N = sizes[0];
  
  
  hsize_t dim2 = sizes[1]; //should be 2
  hid_t posSlabId =  H5Screate_simple(1, &dim2, NULL);


  hid_t QsetId = H5Dopen(groupId,"Q");
  hid_t QspaceId = H5Dget_space(QsetId );