spbound.cpp

pic 模拟程序！面向对象

//--------------------------------------------------------------------
//
// File:	spbound.cpp
//
// Purpose:	Implementation of the SpatialRegionBoundary - a boundary
//		between two spatial regions.  Has methods for passing
//		data between regions.
//
// Version:	$Id: spbound.cpp,v 1.40 2004/10/08 22:32:36 yew Exp $
//
//--------------------------------------------------------------------

/*
====================================================================
Revision History

Cary (23 Jan 00) Added header, implementation of recvParticles, which
	is like emit, except that one does not further weighting of
	the particles.  Did some indentation for readability.
CVS1.27.2.4 (Cary 26 Jan 00) Reorganized to put constructor first.
	Not necessary, I see now, but made sure that j indexing
	was increasing.  Implemented sendStripe and recvStripe.
	Still not checked as we have no way of looking at MPI runs.

(PeterM 1/15/01) Deeper look, allowing Enode/Bnode continuity.
====================================================================
*/

/* These are template defines.  Make it easier not to screw things up.  */

/*The following advance the E field components 1 2 and 3 one timestep.
 "down" means lower index k.  "here" means same index.  "right" means higher index j.
	"left" means lower index j "up" means higher index k*/
#define NEW_E1(iC,iLhere,iLdown,Bhere,Bdown,Ihere) ( (iC).e1() * ( (Bhere).e3()*(iLhere).e3() \
		- (Bdown).e3()*(iLdown).e3() -(Ihere).e1()))


#define NEW_E2(iC,iLhere,iLleft,Bhere,Bleft,Ihere) ( (iC).e2() * ( -(Bhere).e3()*(iLhere).e3() \
	+ (Bleft).e3()*(iLleft).e3() - (Ihere).e2()))

#define NEW_E3(iC,iLhere,iLdown,iLleft,Bhere,Bdown,Bleft,Ihere)  ( (iC).e3() * ( (Bdown).e1()*(iLdown).e1()\
	- (Bhere).e1()*(iLhere).e1()\
	+ (Bhere).e2()*(iLhere).e2()\
	- (Bleft).e2()*(iLleft).e2() - (Ihere).e3()))

/* The following advance the B field components 1 2 and 3 one timestep:  directions are as above. */

#define NEW_B1(Ehere,Eup)        ( (Eup).e3() - (Ehere).e3())
#define NEW_B2(Ehere,Eright)     ( (Ehere).e3() - (Eright).e3())
#define NEW_B3(Ehere,Eright,Eup)  ( (Ehere).e1() - (Eup).e1() + (Eright).e2() - (Ehere).e2())





#include "spbound.h"
#include "fields.h"
#include "mapNGDs.h"


#ifdef _MSC_VER
using std::cout;
using std::cerr;
using std::endl;
#endif

SpatialRegionBoundary::SpatialRegionBoundary(oopicList <LineSegment> *segments) throw (Oops) :
  Boundary(segments) 
   {

	 if(segments->nItems() > 1) {
     stringstream ss (stringstream::in | stringstream::out);
     ss << "SpatialRegionBoundary::SpatialRegionBoundary: Error: \n"<<
      "Spatial region boundaries CANNOT have multiple segments.\n";

		std::string msg;
    ss >> msg;
    Oops oops(msg);
    throw oops;    // exit() SpatialRegionGroup::constructBoundaryParams:  
    /*"SpatialRegionGroup::makeSRB: 
    SpatialRegionBoundaryParams::CreateCounterPart: */

	 }
	 
	 if(j1!=j2 && k1!=k2) {
     stringstream ss (stringstream::in | stringstream::out);
     ss << "SpatialRegionBoundary::SpatialRegionBoundary: Error: \n"<<
		    "Spatial region boundaries MUST be either horizontal or vertical.\n";
		 
     std::string msg;
     ss >> msg;
     Oops oops(msg);
     throw oops;    // exit() SpatialRegionGroup::constructBoundaryParams:
    /*"SpatialRegionGroup::makeSRB: 
    SpatialRegionBoundaryParams::CreateCounterPart: */

	 }

	 int length = MAX(abs(k1-k2),abs(j1-j2));

	 // memory allocation is 1 larger to allow guard cells for a
	 // simpler field solve.
	 lJa = new Vector3[length+1];  // closest guard layer
	 lJb = new Vector3[length+1];  // next closest guard layer
	 lJc = new Vector3[length+1];  // furthest guard layer
	 lEa = new Vector3[length+1];  // closest guard layer
	 lEb = new Vector3[length+1];  // next closest guard layer
	 lEc = new Vector3[length+1];  // furthest guard layer
	 lBa = new Vector3[length+1];  // closest guard layer
	 lBb = new Vector3[length+1];  // next closest guard layer
	 lBc = new Vector3[length+1];  // furthest guard layer
	 liLa = new Vector3[length+1];  // closest guard layer
	 liLb = new Vector3[length+1];  // next closest guard layer
	 liLc = new Vector3[length+1];  // furthest guard layer

	 /* leave the old variables so it will compile */
	 lJ = new Vector3[length+1];  // soon obsolete
	 lE = new Vector3[length+1];  // soon obsolete
	 lB = new Vector3[length+1];  // soon obsolete
	 liL = new Vector3[length+1]; // soon obsolete

	 // Sending 3 layers for simplicity
	 lJsend = new Vector3[3*length+4];
	 lEsend = new Vector3[3*length+4];
	 lBsend = new Vector3[3*length+4];

         //
         // the following buffers are needed for the communication of
         // the neutral gas density (in NGD objects) data strucures
         // across boundaries; memory for them will be allocated only 
         // if there are NGDs; dad 05/16/01
         // 
         numElementsNGDbuf = 0; 
         lNGDbuf = 0;
         lNGDbufSend = 0;
         
	 //zero the memory

	 memset(lJa,0,(length+1) * sizeof(Vector3));
	 memset(lEa,0,(length+1) * sizeof(Vector3));
	 memset(lBa,0,(length+1) * sizeof(Vector3));
	 memset(liLa,0,(length+1) * sizeof(Vector3));

	 memset(lJb,0,(length+1) * sizeof(Vector3));
	 memset(lEb,0,(length+1) * sizeof(Vector3));
	 memset(lBb,0,(length+1) * sizeof(Vector3));
	 memset(liLb,0,(length+1) * sizeof(Vector3));

	 memset(lJc,0,(length+1) * sizeof(Vector3));
	 memset(lEc,0,(length+1) * sizeof(Vector3));
	 memset(lBc,0,(length+1) * sizeof(Vector3));
	 memset(liLc,0,(length+1) * sizeof(Vector3));

	 memset(lJ,0,(length+1) * sizeof(Vector3));
	 memset(lE,0,(length+1) * sizeof(Vector3));
	 memset(lB,0,(length+1) * sizeof(Vector3));
	 memset(liL,0,(length+1) * sizeof(Vector3));

	 memset(lJsend,0,(3*length+4) * sizeof(Vector3));
	 memset(lEsend,0,(3*length+4) * sizeof(Vector3));
	 memset(lBsend,0,(3*length+4) * sizeof(Vector3));
	 theLink = 0;
	 BCType=SPATIAL_REGION_BOUNDARY;
	 ptclPassArraySize=100;
	 ptclPassArray = new particlePassDat[ptclPassArraySize];
	 number_Particles_passing=0;
	 

}

// We need to set up a special iC for the Fields:: so it doesn't interfere
// with the special field solve we do for E2 and E3 on the boundary.

void SpatialRegionBoundary::setPassives() {
  int j,k;
  if(j1==j2) // vertical boundary, set iC2 and iC3 to zero
	 {
		for(k=k1;k<=k2;k++) {
		  fields->set_iL1(j1,k,2*fields->getiL()[j1][k].e1());  // fix iL: default is wrong.
		  fields->set_iC2(j1,k,0);
		  fields->set_iC3(j1,k,0);
		}
	 }
  else  // horizontal boundary, set iC1 and iC3 to zero
	 {
		for(j=j1;j<=j2;j++) {
		  fields->set_iL2(j,k1,2*fields->getiL()[j][k1].e2());  // fix iL: default is wrong.
		  fields->set_iC1(j,k1,0);
		  fields->set_iC3(j,k1,0);
		}
	 }
}
		  

//  This is the function which will assemble and send the currents
//  and fields necessary for the other boundary to solve for its
//  fields.  By convention for this boundary, normal points IN

void SpatialRegionBoundary::putCharge_and_Current(Scalar t) {

  Vector3 **intBdS = fields->getIntBdS();
  Vector3 **intEdl = fields->getIntEdl();
  Vector3 **I = fields->getI();
  // Vector3 **iC = fields->getiC();
  // Vector3 **iL = fields->getiL();

// For indexing
  int kl = MIN(k1,k2);
  int kh = MAX(k1,k2);
  int jl = MIN(j1,j2);
  int jh = MAX(j1,j2);
  int length;

// Begin the request
  if(!theLink) return;
  theLink->askFields();
  length = MAX(abs(k2-k1),abs(j2-j1));

// Horizontal boundary data send
  if(k1==k2) for(int j=jl;j<jh;j++) {  

    if(normal==-1) {
// Real region below
      lBsend[j-jl]=Vector3(intBdS[j][k1-1].e1(),0,intBdS[j][k1-1].e3());
      lEsend[j-jl]=intEdl[j][k1-1];
      lJsend[j-jl]=Vector3(I[j][k1].e1(),I[j][k1-1].e2(),I[j][k1].e3());
    }

    else {
// Real region above
      lBsend[j-jl]=Vector3(intBdS[j][k1].e1(), 0, intBdS[j][k1].e3());
      lEsend[j-jl]=Vector3(intEdl[j][k1].e1(), intEdl[j][k1+1].e2(),
		intEdl[j][k1+1].e3());
      lJsend[j-jl]=Vector3(I[j][k1].e1(), I[j][k1].e2(), I[j][k1].e3());
    }

  }

// Vertical boundary data send
  else for(int k=kl;k<=kh;k++) {
	  int kk = k - kl;

// Real region left
    if(normal==-1) {
/*      lBsend[k-kl]=Vector3(0,intBdS[j1-1][k].e2(),intBdS[j1-1][k].e3());
      lEsend[k-kl]=intEdl[j1-1][k];
      lJsend[k-kl]=Vector3(I[j1-1][k].e1(),I[j1][k].e2(),I[j1][k].e3()); */
		 //New way of doing this, send 3 full layers
		 // closest layer first
		 lBsend[kk] = intBdS[j1][k];
		 lBsend[length + 1 + kk] = intBdS[j1-1][k];
		 lBsend[2*length +2 + kk] = intBdS[j1-2][k];

		 lEsend[kk] = intEdl[j1][k];
		 lEsend[length + 1 + kk] = intEdl[j1-1][k];
		 lEsend[2*length +2 + kk] = intEdl[j1-2][k];

		 lJsend[k-kl] = I[j1][k];
		 lJsend[length + 1 + kk] = I[j1-1][k];
		 lJsend[2*length +2 + kk] = I[j1-2][k];
		 
    }

// Real region right
    else {
      /*lBsend[k-kl]=Vector3(0,intBdS[j1][k].e2(),intBdS[j1][k].e3());
      lEsend[k-kl]=Vector3(intEdl[j1][k].e1(), intEdl[j1+1][k].e2(),
		intEdl[j1+1][k].e3());
      lJsend[k-kl]=I[j1][k]; */
		 lBsend[kk] = intBdS[j1][k];
		 lBsend[length + 1 + kk] = intBdS[j1+1][k];
		 lBsend[2*length +2 + kk] = intBdS[j1+2][k];

		 lEsend[kk] = intEdl[j1][k];
		 lEsend[length + 1 + kk] = intEdl[j1+1][k];
		 lEsend[2*length +2 + kk] = intEdl[j1+2][k];

		 lJsend[kk] = I[j1][k];
		 lJsend[length + 1 + kk] = I[j1+1][k];
		 lJsend[2*length +2 + kk] = I[j1+2][k];

    }

  }

/*
#ifdef MPI_VERSION
  cout << "putCharge_and_Current sending fields:" << endl;
  cout << "E1 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lEsend[kk].e1();
  cout << endl;
  cout << "E2 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lEsend[kk].e2();
  cout << endl;
  cout << "E3 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lEsend[kk].e3();
  cout << endl;
  cout << "B1 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lBsend[kk].e1();
  cout << endl;
  cout << "B2 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lBsend[kk].e2();
  cout << endl;
  cout << "B3 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lBsend[kk].e3();
  cout << endl;
#endif
*/

// Send the fields through the link
  theLink->sendFields(lEsend,lBsend,lJsend);
}                  

// This is the function which does the fieldsolve.
// It does it in such a way that toNodes will be wrong,
// and it does it in such a way that serial operation is
// possible.
// By convention for this boundary, normal points IN

void SpatialRegionBoundary::applyFields(Scalar t,Scalar dt) {

  Scalar dt2 = dt/2;
  Vector3 **intBdS = fields->getIntBdS();
  Vector3 **intEdl = fields->getIntEdl();
  Vector3 **I = fields->getI();
  Vector3 **iC = fields->getiC();
  Vector3 **iL = fields->getiL();
 
  int kl = MIN(k1,k2);
  int kh = MAX(k1,k2);
  int jl = MIN(j1,j2);
  int jh = MAX(j1,j2);
  int K = fields->getK();
  int J = fields->getJ();
  if(!theLink) return;
  // old version
  //theLink->waitFields();  //  do a synchronize
  theLink->waitFields(lEa,lEb,lEc,lBa,lBb,lBc,lJa,lJb,lJc);

  // Advance B 1/2 step

  if(dt > 0) {
     SRBadvanceB(t,dt);

     // Advance E 1 step
     SRBadvanceE(t,dt);

  
     // Advance B 1/2 step is now NECESSARY for weighting to the nodes properly.
     
     SRBadvanceB(t,dt);
  }

}



//  This function sends the iL's the paired boundary will need to
//  complete the field solve.
void SpatialRegionBoundary::send_liL() {
  int j,k;
  Vector3 **iL = fields->getiL();
  if(k1==k2) {
	 for(j=j1;j<j2;j++) {
		if(normal==1)  //send the iL's from above
		  lBsend[j-j1] = iL[j][k1+1];
		else //send the iL's from below
		  lBsend[j-j1] = iL[j][k1-1];
	 }
  }
  else {
	 for(k=k1;k<k2;k++) {
		if(normal==1)  //send the iL's from right
		  lBsend[k-k1] = iL[j1+1][k];
		else //send the iL's from left
		  lBsend[k-k1] = iL[j1-1][k];
	 }
  }
  //  Actually send the data over.
  theLink->askFields();  //initialize the transfer
  theLink->sendFields(lEsend,lBsend,lJsend); // do our send
  theLink->waitFields();  //synchronize
}

//  This function initializes the passives across the SRB's.
void SpatialRegionBoundary::initialize_passives() {
  int j,k;
  Vector3 **iL = fields->getiL();
  if(k1==k2) {
	 for(j=j1;j<=j2;j++) {
		if(normal==1)  //send the iL's from below
		  lBsend[j-j1] = iL[j][k1+1];
		else //send the iL's from above
		  lBsend[j-j1] = iL[j][k1-1];
	 }
  }
  else {
	 for(k=k1;k<=k2;k++) {
		if(normal==1)  //send the iL's from right
		  lBsend[k-k1] = iL[j1+1][k];
		else //send the iL's from left
		  lBsend[k-k1] = iL[j1-1][k];
	 }
  }
  //  Actually send the data over.
  theLink->askFields();  //initialize the transfer
  theLink->sendFields(lEsend,lBsend,lJsend); // do our send
  theLink->waitFields();  //synchronize
  get_liL();
}

//  This function sends the iL's the paired boundary will need to
//  complete the field solve.
void SpatialRegionBoundary::get_liL() {
  int j,k;
    if(k1==k2) {
	 for(j=j1;j<=j2;j++) {
		if(normal==1) {
		  liL[j-j1] = lB[j-j1];
		  lB[j-j1]=0;
		}
		else {
		  liL[j-j1] = lB[j-j1];
		  lB[j-j1]=0;
		}
	 }
  }
  else {
	 for(k=k1;k<k2+1;k++) {
		if(normal==1) { 
		  liL[k-k1] = lB[k-k1];
		  lB[k-k1]=0;
		}
		else {
		  liL[k-k1] = lB[k-k1];
		  lB[k-k1]=0;
		}
	 }
  }
}

void SpatialRegionBoundary::toNodes()
{
  oopicListIter <LineSegment> lsI(*segments);
  for(lsI.restart();!lsI.Done();lsI++) {
     int j1,j2,k1,k2,normal,*points;  //local copies of above
     j1=(int)lsI()->A.e1();
     k1=(int)lsI()->A.e2();
     j2=(int)lsI()->B.e1();
     k2=(int)lsI()->B.e2();
     normal = lsI()->normal;
     points=lsI()->points; 
	 
     int j, k,K;
     Scalar w1m=0.5, w1p=0.5;
     Vector3** ENode = fields->getENode();
     Vector3** BNode = fields->getBNode();
     Vector3** intEdl = fields->getIntEdl();
     Vector3** intBdS = fields->getIntBdS();
     Grid* grid = fields->get_grid();
     K = grid->getK();
     BNode[j1][0] = 0;  // unbounded edges

     if(j1==j2 && k1==k2) return;  // don't do toNodes on one-point boundaries.
     if (j1==j2)            // vertical boundary
	{
	   for (j=j1, k=k1; k<=k2; k++) {
	      int kk = k - k1;
	      if(normal == 1) {
		 // ASSUME:  Uniform mesh

		 Scalar w2m = grid->dl2(j,k) / (grid->dl2(j,k) + grid->dl2(j,MAX(0,k-1)));
		 Scalar w2p = 1-w2m;
		 Scalar w1m = 0.5;  // uniform mesh assumption
		 Scalar w1p = 0.5;  // uniform mesh assumption