dielectr.cpp

pic 模拟程序！面向对象

// dump boundary type as attribute
  dims = 1;
  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

 rank = 2;
  hsize_t     dimsf2[2];
  //  Scalar data[nPoints][2];
  dimsf2[0] = nPoints;
  dimsf2[1] = 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 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);}
      x = X[j][k].e1();
      y = X[j][k].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);


  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 DielectricRegion::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");

   Scalar x,y;
   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();

   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);}
      x = X[j][k].e1();
      y = X[j][k].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 DielectricRegion::restoreH5(dumpHDF5 &restoreObj,int bType,string whichBoundary,int nQuads)
{
  //cerr << "entering DielectricRegion::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());
    


  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.
  
// Open the dataset for reading.  Determine shape and set field size.





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

  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);
    
    int kmax = MAX(k1,k2);
    int jmax = MAX(j1,j2);
    int jmin = MIN(j1,j2);
    int kmin = MIN(k1,k2);
    Q[getIndex(jl,kl)]+= lQ;
  }	
   
  H5Sclose(dataspaceId);
  H5Dclose(datasetId);
  H5Sclose(QspaceId);
  H5Dclose(QsetId);
  H5Sclose(posSlabId);
  
  
  H5Gclose(groupId);
  H5Fclose(fileId);

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

int DielectricRegion::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);

      int kmax = MAX(k1,k2);
      int jmax = MAX(j1,j2);
      int jmin = MIN(j1,j2);
      int kmin = MIN(k1,k2);
      Q[getIndex(jl,kl)]+= lQ;

      //find which edge it is on and addQ it
      /*		if(jl == jmin) {
			if(kl <=kmax && kl >=kmin) {
                        //addQ(LEFT,jl,kl,lQ);

			}
                        }

                        if(jl==jmax) {
			if(kl <=kmax && kl >=kmin) {
                        addQ(RIGHT,jl,kl,lQ);
                        }
                        }

                        if(kl == kmin) {
			if(jl <=jmax && jl >=jmin) {
                        addQ(BOTTOM,jl,kl,lQ);
			}
                        }
                        if(kl == kmax) {
			if(jl <=jmax && jl >=jmin) {
                        addQ(TOP,jl,kl,lQ);
			}
                        }

                        */
   }	
   //#endif
   return(TRUE);


}

int DielectricRegion::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);
}

DielectricTriangle::
DielectricTriangle(oopicList <LineSegment> *segments,
                   Scalar _er, int _QuseFlag) :Dielectric(segments,_er,0)
{
   QuseFlag=_QuseFlag;
   if(segments->nItems() > 1) {
      cout << "Warning, Dielectric Triangles can only have 1 segment.\n";
      cout << "Check your input file.\n";
   }

   /*	if(j1>j2) {  //swap the two points.
        int jt,kt;
        jt=j2;j2=j1;j1=jt;
        kt=k2;k2=k1;k1=kt;
	}
        */
   if(QuseFlag)
      {
         fprintf(stderr,"DielectricTriangles don't support QuseFlag\n");
         QuseFlag = 0;
      }
   epsilon = _er/iEPS0;

   BoundaryType = normal*(DIELECTRIC_TRIANGLE_UP);
   BCType = DIELECTRIC_BOUNDARY;
}

void DielectricTriangle::setPassives()
{
   /*
      int j,k,jl,kl,jh,kh;
      int J=fields->getJ();
      int K=fields->getK();
      int *points = segments->head->data->points;
      // If the normal is positive, fill the triangle upwards
      kh=MAX(k1,k2);
      if(normal > 0) {
      for(j=0;j<4*abs(j2-j1)+4;j+=2) {
      kl = points[j+1];
      jl = points[j];
      for(k=kl;k<=kh;k++)
      {
      if (jl==0 || jl==J || k==0 || k==K)
      fields->getPoissonSolve()->set_coefficient(jl,k,NEUMANN,fields->get_grid());
      else
      fields->getPoissonSolve()->set_coefficient(jl,k,FREESPACE,fields->get_grid());
      }
      }
      }
      else {
      kl = MIN(k1,k2);
      for(j=0;j<4*abs(j2-j1)+4;j+=2) {
      kh = points[j+1];
      jl = points[j];
      for(k=kl;k<=kh;k++)
      {
      if (jl==0 || jl==J || k==0 || k==K)
      fields->getPoissonSolve()->set_coefficient(jl,k,NEUMANN,fields->get_grid());
      else
      fields->getPoissonSolve()->set_coefficient(jl,k,FREESPACE,fields->get_grid());
      }
      }
      }
      */
}

void DielectricTriangle::collect(Particle &p, Vector3& dxMKS)
{
   delete &p;
}

//set the permittivity inside the triangle

void DielectricTriangle::setEpsilon()
{
   int j,k,jl,kl,kh;

   int *points = segments->head->data->points;
   // If the normal is positive, fill the triangle upwards
   kh=MAX(k1,k2);
   if(normal > 0) {
      for(j=0;j<4*abs(j2-j1)+4;j+=2) {
         kl = points[j+1];
         jl = points[j];
         for(k=kl;k<=kh;k++)
            fields->set_epsi(epsilon,jl,k);
      }
   }
   else {
      kl = MIN(k1,k2);
      for(j=0;j<4*abs(j2-j1)+4;j+=2) {
         kh = points[j+1];
         jl = points[j];
         for(k=kl;k<=kh;k++)
            fields->set_epsi(epsilon,jl,k);
      }
   }
}	

void DielectricTriangle::setBoundaryMask(Grid &grid) {

   Boundary::setBoundaryMask(grid);  //this does a lot of the work we need to do.

   //for the other two boundaries of the triangle, we need to figure out where to
   // put them.
   int jm,km,jh,kh,jl,kl;
   kh = MAX(k1,k2);
   jh = MAX(j1,j2);
   jl = MIN(j1,j2);
   kl = MIN(k1,k2);
   
   if(normal > 0) {
      km = kh;
      jm = ( k1 < k2) ? jl:jh;
   }
   else {
      km=kl;
      jm = (k1 < k2)? jh:jl;
   }

   //Horizontal boundary
   for(int j=jl;j<jh;j++) grid.setBoundaryMask1(j,km,this,0);
   
   //Vertical boundary
   for(int k=kl;k<kh;k++) grid.setBoundaryMask2(jm,k,this,0);
}




int DielectricTriangle::Dump(FILE * DMPFile)
{
   return(TRUE);
}