dielectr.cpp

pic 模拟程序！面向对象

  for(i=0;i<nQuads;i++) {	
    Scalar x,y,lQ; //dummy;
    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];
      
    offset[0]  = 0;  count[0]  = 1;
    status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, offset, NULL, count, NULL);
    offset[0]  = i; offset[1]  = 0;
    count[0]  = 1; count[1]  = 1;     
    status = H5Sselect_hyperslab(QspaceId, H5S_SELECT_SET, offset, NULL, count, NULL);
    status = H5Dread(QsetId, scalarType, posSlabId, QspaceId, H5P_DEFAULT, &lQ);
    

// 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)
               Q[ki]= lQ;

         }
      else   // horizontal
         {
            int jmax = abs(j2-j1) + 1;
            int ji = jl - MIN(j1,j2);
            if(ji <= jmax && ji >=0)
               Q[ji]= lQ;
         }
   }	


  H5Sclose(dataspaceId);
  H5Dclose(datasetId);
  H5Sclose(QspaceId);
  H5Dclose(QsetId);
  H5Sclose(posSlabId);
  
  H5Gclose(groupId);
  H5Fclose(fileId);

  delete [] sizes;  
  return(1);
}
#endif


int Dielectric::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;
   
   // OK, this dielectric must be us.
   

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

   for(i=0;i<nQuads;i++) {	
      Scalar x,y,lQ,dummy;
      Vector2 GridLoc;
      XGRead(&x,ScalarInt,1,DMPFile,ScalarChar);
      XGRead(&y,ScalarInt,1,DMPFile,ScalarChar);
      XGRead(&lQ,ScalarInt,1,DMPFile,ScalarChar);
      XGRead(&dummy,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)
               Q[ki]= lQ;

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


int Dielectric::Restore_2_00(FILE * DMPFile, int j1test,
                             int k1test, int j2test, int k2test)
{
   if ((j1 == j1test)&&(k1 == k1test)&&(j2 == j2test)&&(k2 == k2test))

      {
         //#ifdef UNIX
         for(int i=0; i<nPoints; i++)
            XGRead(&Q[i], ScalarInt, 1, DMPFile, ScalarChar);
         //#endif
         return(TRUE);

      }
   else
      return(FALSE);
}

//===============================================================
// DielectricRegion.  Must construct Q with different array size than
//	Dielectric, so pass QuseFlag=0 to prevent allocation in Dielectric()

DielectricRegion::DielectricRegion(oopicList <LineSegment> *segments, 
                                   Scalar er, int _QuseFlag, Scalar reflection, Scalar refMaxE) 
: Dielectric(segments, er, 0)
{

   if(segments->nItems() > 1) {
      cout << "Warning, Dielectric Regions can only have 1 segment.\n";
      cout << "Check your input file.\n";
   }

   QuseFlag = _QuseFlag;					// must set here since we send 0 to Dielectric()
   if (QuseFlag)
      {
         nPoints = 2*(j2 - j1 + k2 - k1);
         Q = new Scalar[nPoints*2];
         memset(Q, 0, nPoints*sizeof(Scalar));
      }
   BoundaryType = DIELECTRIC_REGION;
   BCType = DIELECTRIC_BOUNDARY;
}

DielectricRegion::~DielectricRegion()
{
}

void DielectricRegion::setPassives()
{
}

#define DELTA 1e-6

void DielectricRegion::collect(Particle& p, Vector3& dxMKS)
{
   int j=0;
   EDGE _edge;
   Scalar xIntercept=0.;
   // Compute the edge which collected:
   if (fabs(p.get_x().e1() - j1) < DELTA)	//	left edge
      {
         xIntercept = p.get_x().e2();
         _edge = LEFT;
         _delta.set_e1(j1 ? -DELTA*j1 : -DELTA);
         _delta.set_e2(0);
         _unitNormal.set_e1(-1);
         _unitNormal.set_e2(0);
      }
   else if (fabs(p.get_x().e1() - j2) < DELTA)	// right edge
      {
         xIntercept = p.get_x().e2();
         _edge = RIGHT;
         _delta.set_e1(DELTA*j2);
         _delta.set_e2(0);
         _unitNormal.set_e1(1);
         _unitNormal.set_e2(0);
      }
   else if (fabs(p.get_x().e2() - k1) < DELTA)	// bottom edge
      {
         xIntercept = p.get_x().e1();
         _edge = BOTTOM;
         _delta.set_e1(0);
         _delta.set_e2(k1 ? -DELTA*k1 : DELTA);
         _unitNormal.set_e1(0);
         _unitNormal.set_e2(-1);
      }
   else if (fabs(p.get_x().e2() - k2) < DELTA)	// top edge
      {
         xIntercept = p.get_x().e1();
         _edge = TOP;
         _delta.set_e1(0);
         _delta.set_e2(DELTA*k2);
         _unitNormal.set_e1(0);
         _unitNormal.set_e2(1);
      }
   else {
      cout << "DielectricRegion::collect failed to find an edge" << endl;
      printf("delta=(%g,%g)   x=(%g,%g)   s=%s   j=%d, k=%d\n", 
	_delta.e1(), _delta.e2(), p.get_x().e1(), p.get_x().e2(), 
	p.get_speciesPtr()->get_name().c_str(), j1, k1);
// JRC: if no edge, cannot continue.
      return;
   }
   if (QuseFlag)
      {
         j = (int) xIntercept;
         addQ(_edge, j, xIntercept - j, p.get_q());
      }
   Boundary::collect(p, dxMKS); // for stats
   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 && n>0){
            Scalar q = -n*p.get_np2c()*sIter()->get_q();
            addQ(_edge, j, xIntercept - j, q);
         }
      }
   }
   delete &p;
}

//-------------------------------------------------------------------
// Set the permittivity inside the dielectric region
void DielectricRegion::setEpsilon()
{
   //return;
   for (int j = j1; j < j2; j++)
      for (int k = k1; k < k2; k++)
         fields->set_epsi(epsilon, j, k);
}

//-------------------------------------------------------------------
// Add q to the surface charge array; must be mapped from perimeter
// to linear array.
void DielectricRegion::addQ(EDGE edge, int i, Scalar w, Scalar q)
{
   int index1=0;
   int index2=0;
   switch (edge)
      {
       case BOTTOM:
          index1 = i - j1;
          index2 = index1 + 1;
          break;
        case RIGHT:
           index1 = i - k1 + j2 - j1;
          index2 = index1 + 1;
          break;
        case TOP:
           index1 = 2*j2 - j1 + k2 - k1 - i;
          index2 = index1 - 1;
          break;
        case LEFT:
           index1 = 2*(j2 - j1) + 2*k2 - k1 - i;
          index2 = index1 - 1;
          index1 %= nPoints;				// wrap around to 0
          break;             
       }
   Q[index1] += q*(1 - w);
   Q[index2] += q*w;
}

int DielectricRegion::getIndex(int j, int k)
{
   if (k == k1) return j - j1;    				//	bottom
   else if (j == j2)	return k - k1 + j2 - j1;		//	right
   else if (k == k2) return 2*j2 - j1 + k2 - k1 - j;	//	top
   else if (j == j1) return 2*(j2 - j1 + k2) - k1 - k;	//	left
   else exit (-1);							// points must lie on surface
   return j-j1;  //  goofy default case.
}

void DielectricRegion::putCharge_and_Current(Scalar t)
{
   if (!QuseFlag) return;
   Scalar **rho=fields->getrho();
   Scalar **cellvolume = fields->get_grid()->get_halfCellVolumes();
   for (int i=0; i < nPoints; i++)
      {
         int j, k;
         if (i < j2 - j1) {j = j1 + i; k = k1;}
         else if (i < j2 - j1 + k2 - k1) {j = j2; k = i - j2 + j1 + k1;}
         else if (i < 2*(j2 - j1) + k2 - k1) {j = 2*j2 - i - j1 + k2 - k1; k = k2;}
         else {j = j1; k = 2*k2 - k1 - i + 2*(j2 - j1);}
         rho[j][k] += Q[i]/cellvolume[j][k];
      }
}

void DielectricRegion::putCharge(void)
{
   if (!QuseFlag) return;
   Scalar **Charge=fields->getQ();
   for (int i=0; i < nPoints; i++)
      {
         int j, k;
         if (i < j2 - j1) {j = j1 + i; k = k1;}
         else if (i < j2 - j1 + k2 - k1) {j = j2; k = i - j2 + j1 + k1;}
         else if (i < 2*(j2 - j1) + k2 - k1) {j = 2*j2 - i - j1 + k2 - k1; k = k2;}
         else {j = j1; k = 2*k2 - k1 - i + 2*(j2 - j1);}
         Charge[j][k] += Q[i];
      }
}

//---------------------------------------------------------------------
//	setBoundaryMask sets the boundary masks in Grid for determing boundaries.
// Boundary::setBoundaryMask() is overloaded to deal with volumetric 
// boundaries.

void DielectricRegion::setBoundaryMask(Grid &grid)
{
   int j, k=k2;
   // setBoundaryMask might need to be replaced with SetNodeBoundary so
   // particles can live in a dielectric. 
   // ie. M. Turner changing epsilon speed up.
   /*
      for (j=j1; j<j2; j++)
      for (int k=k1; k<=k2; k++)
      grid.setBoundaryMask1(j, k, this);

      for (j=j1; j<=j2; j++)
      for (int k=k1; k<k2; k++)
      grid.setBoundaryMask2(j, k, this);
      */
   for (j=j1; j<j2; j++)
      for (k=k1; k<k2; k++){
         grid.setCellMask(j, k, this);
         grid.setBoundaryMask1(j, k, this, 0);
         grid.setBoundaryMask2(j, k, this, 0);
      }

   for (j=j1; j<j2; j++)
      grid.setBoundaryMask1(j, k, this, 0);
   
   for (k=k1; k<k2; k++)
      grid.setBoundaryMask2(j, k, this, 0);

}

//---------------------------------------------------------------------
//	Weight the fields to nodal locations along the edges of the dielectric
//	region.  Must be handled slightly differently than Boundary::toNodes()
//	just to account for four edges.  Assume only normal fields need be
//	updated.

void DielectricRegion::toNodes()
{
   Vector3** intEdl = fields->getIntEdl();
   Vector3** ENode = fields->getENode();
   Grid* grid = fields->get_grid();
   for (int j = j1; j <= j2; j++)		// along x1
      {
         if (k1 > 1) ENode[j][k1].set_e2(intEdl[j][k1-1].e2()/grid->dl2(j, k1-1));
         if (k2 < grid->getK())	ENode[j][k2].set_e2(intEdl[j][k2].e2()/grid->dl2(j, k2));
         /*		if (k1 > 0)	ENode[j][k1].set_e2((epsilon*intEdl[j][k1].e2()/grid->dl2(j, k1) 
			- Q[getIndex(j, k1)]/grid->dS2(j,k1))/fields->get_epsi(j, k1-1));
                        if (k2 < grid->getK())	ENode[j][k2].set_e2((Q[getIndex(j, k2)]/grid->dS2(j,k1)
			- epsilon*intEdl[j][k2-1].e2()/grid->dl2(j, k2-1))/fields->get_epsi(j, k2));
                        */	}        
   for (int k = k1; k <= k2; k++)		// along x2
      {
         if (j1 > 0) ENode[j1][k].set_e1(intEdl[j1-1][k].e1()/grid->dl1(j1-1, k));
         if (j2 < grid->getJ()) ENode[j2][k].set_e1(intEdl[j2][k].e1()/grid->dl1(j2, k));
         /*		if (j1 > 0) ENode[j1][k].set_e1((epsilon*intEdl[j1][k].e1()/grid->dl1(j1, k)
			- Q[getIndex(j1, k)]/grid->dS1(j1, k))/fields->get_epsi(j1-1, k));
                        if (j2 < grid->getJ()) ENode[j2][k].set_e1((Q[getIndex(j2, k)]/grid->dS1(j2, k)
			- epsilon*intEdl[j2-1][k].e1()/grid->dl1(j2-1, k))/fields->get_epsi(j2, k));
                        */	}
}
#ifdef HAVE_HDF5
int DielectricRegion::dumpH5(dumpHDF5 &dumpObj,int number)
{

  //cerr << "Entered DielectricRegion::Dump--HDF5.\n";
  hid_t fileId, groupId, datasetId,dataspaceId;
  hid_t  attrdataspaceId,attributeId;
  herr_t status;
  string outFile;
  hsize_t rank; 
  hsize_t     dimsf[1];   
  hsize_t dims;
  hid_t scalarType = dumpObj.getHDF5ScalarType();
#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 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=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();
  
  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);