load.cpp

pic 模拟程序！面向对象

						gridded_density[j][k] =0;
				}
			cout << "finished reading density file " << den_file << endl;
			fclose (openfile);
			
		}
	else
		{
      stringstream ss (stringstream::in | stringstream::out);
      ss << "Load::init_file_load: Error: \n"<<
        "Density file not valid." << endl;
			std::string msg;
      ss >> msg;
      Oops oops(msg);
      throw oops;    // exit()  Load::Load

		}
}

void Load::IntegrateGriddedLoad(void)  //and set back ground rho and show initial load
{

	int j,k;
	Scalar total_density=0;
	Scalar max_density=0;
	Scalar ** int_density;
	int_density = new Scalar *[J+1];
	for (j=0; j<=J; j++){
		int_density[j] = new Scalar[K + 1];
		for(k=0;k<=K;k++){ 
			int_density[j][k]=0;
		}
	}

	Vector2 p1Grid = grid->getGridCoords(p1);
	Vector2 p2Grid = grid->getGridCoords(p2)+Vector2(1,1);

	/***doing a second order integration of the gridded_density***/

	Boundary*** CellMask = grid->GetCellMask();
	Scalar dV=0;
	Scalar fx, fy;
	Scalar fx13,fx31,fx12;
	Scalar fy13,fy31,fy12;

	for (j=(int)p1Grid.e1(); j<MIN((int)p2Grid.e1(),J); j++)
		for (k=(int)p1Grid.e2(); k<MIN((int)p2Grid.e2(),K); k++){ 	
			if (CellMask[j][k]==NULL){ //FREESPACE
				// first if is getting all the half filled cells
				if (((j==(int)p1Grid.e1())&&(p1Grid.e1()!=0))||
						((j==((int)p2Grid.e1()-1))&&((p1Grid.e1()-1)!=J))||
						((k==((int)p1Grid.e2())&&(p1Grid.e2()!=0)))||
						((k==((int)p2Grid.e2()-1))&&((p2Grid.e2()-1)!=K))){  //this is over checking, trying to increase speed
					if ((j==(int)p1Grid.e1())&&(k==(int)p1Grid.e2())){  //lower left corner
						dV = grid->cellVolume(j,k);
						fx = p1Grid.e1()-(int)p1Grid.e1();
						fy = p1Grid.e2()-(int)p1Grid.e2();
						fx13 = cube(1-fx);
						fx31 = 1-cube(fx);
						fx12 = (1+2*fx)*sqr(1-fx);

						fy13 = cube(1-fy);
						fy31 = 1-cube(fy);
						fy12 = (1+2*fy)*sqr(1-fy);
						
						int_density[j][k] += (4*fx13*fy13*gridded_density[j][k]+
																	2*fx12*fy13*gridded_density[j+1][k]+
																	2*fx13*fy12*gridded_density[j][k+1]+
																	fx12*fy12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fx31*fy13*gridded_density[j+1][k]+
																		2*fx12*fy13*gridded_density[j][k]+
																		2*fx31*fy12*gridded_density[j+1][k+1]+
																		fx12*fy12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fx31*fy31*gridded_density[j+1][k+1]+
																			2*fx31*fy12*gridded_density[j+1][k]+
																			2*fx12*fy31*gridded_density[j][k+1]+
																			fx12*fy12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fx13*fy31*gridded_density[j][k+1]+
																		2*fx12*fy31*gridded_density[j+1][k+1]+
																		2*fx13*fy12*gridded_density[j][k]+
																		fx12*fy12*gridded_density[j+1][k])/36*dV;	

					}
					else if (j==((int)p1Grid.e1())&&(k==((int)p2Grid.e2()-1))){  //upper left corner	
						dV = grid->cellVolume(j,k);
						fx = p1Grid.e1()-(int)p1Grid.e1();
						fy = (int)p2Grid.e2()-p2Grid.e2()+1;
						fx13 = cube(1-fx);
						fx31 = 1-cube(fx);
						fx12 = (1+2*fx)*sqr(1-fx);

						fy13 = cube(1-fy);
						fy31 = 1-cube(fy);
						fy12 = (1+2*fy)*sqr(1-fy);
						
						int_density[j][k] += (4*fx13*fy31*gridded_density[j][k]+
																	2*fx12*fy31*gridded_density[j+1][k]+
																	2*fx13*fy12*gridded_density[j][k+1]+
																	fx12*fy12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fx31*fy31*gridded_density[j+1][k]+
																		2*fx12*fy31*gridded_density[j][k]+
																		2*fx31*fy12*gridded_density[j+1][k+1]+
																		fx12*fy12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fx31*fy13*gridded_density[j+1][k+1]+
																			2*fx31*fy12*gridded_density[j+1][k]+
																			2*fx12*fy13*gridded_density[j][k+1]+
																			fx12*fy12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fx13*fy13*gridded_density[j][k+1]+
																		2*fx12*fy13*gridded_density[j+1][k+1]+
																		2*fx13*fy12*gridded_density[j][k]+
																		fx12*fy12*gridded_density[j+1][k])/36*dV;	


					}
					else if ((j==((int)p2Grid.e1()-1))&&(k==((int)p1Grid.e2()))){  //lower right corner
						dV = grid->cellVolume(j,k);
						fx = (int)p2Grid.e1()-p2Grid.e1()+1;
						fy = p1Grid.e2()-(int)p1Grid.e2();

						fx13 = cube(1-fx);
						fx31 = 1-cube(fx);
						fx12 = (1+2*fx)*sqr(1-fx);

						fy13 = cube(1-fy);
						fy31 = 1-cube(fy);
						fy12 = (1+2*fy)*sqr(1-fy);
						
						int_density[j][k] += (4*fx31*fy13*gridded_density[j][k]+
																	2*fx12*fy13*gridded_density[j+1][k]+
																	2*fx31*fy12*gridded_density[j][k+1]+
																	fx12*fy12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fx13*fy13*gridded_density[j+1][k]+
																		2*fx12*fy13*gridded_density[j][k]+
																		2*fx13*fy12*gridded_density[j+1][k+1]+
																		fx12*fy12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fx13*fy31*gridded_density[j+1][k+1]+
																			2*fx13*fy12*gridded_density[j+1][k]+
																			2*fx12*fy31*gridded_density[j][k+1]+
																			fx12*fy12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fx31*fy31*gridded_density[j][k+1]+
																		2*fx12*fy31*gridded_density[j+1][k+1]+
																		2*fx31*fy12*gridded_density[j][k]+
																		fx12*fy12*gridded_density[j+1][k])/36*dV;	
						
					}
					else if ((j==((int)p2Grid.e1()-1))&&(k==((int)p2Grid.e2()-1))){  //upper right corner
						dV = grid->cellVolume(j,k);
						fx = (int)p2Grid.e1()-p2Grid.e1()+1;
						fy = (int)p2Grid.e2()-p2Grid.e2()+1;

						fx13 = cube(1-fx);
						fx31 = 1-cube(fx);
						fx12 = (1+2*fx)*sqr(1-fx);

						fy13 = cube(1-fy);
						fy31 = 1-cube(fy);
						fy12 = (1+2*fy)*sqr(1-fy);
						
						int_density[j][k] += (4*fx31*fy31*gridded_density[j][k]+
																	2*fx12*fy31*gridded_density[j+1][k]+
																	2*fx31*fy12*gridded_density[j][k+1]+
																	fx12*fy12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fx13*fy31*gridded_density[j+1][k]+
																		2*fx12*fy31*gridded_density[j][k]+
																		2*fx13*fy12*gridded_density[j+1][k+1]+
																		fx12*fy12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fx13*fy13*gridded_density[j+1][k+1]+
																			2*fx13*fy12*gridded_density[j+1][k]+
																			2*fx12*fy13*gridded_density[j][k+1]+
																			fx12*fy12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fx31*fy13*gridded_density[j][k+1]+
																		2*fx12*fy13*gridded_density[j+1][k+1]+
																		2*fx31*fy12*gridded_density[j][k]+
																		fx12*fy12*gridded_density[j+1][k])/36*dV;	

					}
					else if (k==(int)p1Grid.e2()){ //lower edge, no corners
						dV = grid->cellVolume(j,k);
						fy = p1Grid.e2()-(int)p1Grid.e2();

						fy13 = cube(1-fy);
						fy31 = 1-cube(fy);
						fy12 = (1+2*fy)*sqr(1-fy);
						
						int_density[j][k] += (4*fy13*gridded_density[j][k]+
																	2*fy13*gridded_density[j+1][k]+
																	2*fy12*gridded_density[j][k+1]+
																	fy12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fy13*gridded_density[j+1][k]+
																		2*fy13*gridded_density[j][k]+
																		2*fy12*gridded_density[j+1][k+1]+
																		fy12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fy31*gridded_density[j+1][k+1]+
																			2*fy12*gridded_density[j+1][k]+
																			2*fy31*gridded_density[j][k+1]+
																			fy12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fy31*gridded_density[j][k+1]+
																		2*fy31*gridded_density[j+1][k+1]+
																		2*fy12*gridded_density[j][k]+
																		fy12*gridded_density[j+1][k])/36*dV;
						
					}
					else if (j==(int)p1Grid.e1()){  //left edge, no corner
						dV = grid->cellVolume(j,k);
						fx = p1Grid.e1()-(int)p1Grid.e1();
						
						fx13 = cube(1-fx);
						fx31 = 1-cube(fx);
						fx12 = (1+2*fx)*sqr(1-fx);

						int_density[j][k] += (4*fx13*gridded_density[j][k]+
																	2*fx12*gridded_density[j+1][k]+
																	2*fx13*gridded_density[j][k+1]+
																	fx12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fx31*gridded_density[j+1][k]+
																		2*fx12*gridded_density[j][k]+
																		2*fx31*gridded_density[j+1][k+1]+
																		fx12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fx31*gridded_density[j+1][k+1]+
																			2*fx31*gridded_density[j+1][k]+
																			2*fx12*gridded_density[j][k+1]+
																			fx12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fx13*gridded_density[j][k+1]+
																		2*fx12*gridded_density[j+1][k+1]+
																		2*fx13*gridded_density[j][k]+
																		fx12*gridded_density[j+1][k])/36*dV;
						
					}
					else if (k==((int)p2Grid.e2()-1)){  //upper edge, no corner
						dV = grid->cellVolume(j,k);
						fy = (int)p2Grid.e2()-p2Grid.e2()+1;

						fy13 = cube(1-fy);
						fy31 = 1-cube(fy);
						fy12 = (1+2*fy)*sqr(1-fy);
						
						int_density[j][k] += (4*fy31*gridded_density[j][k]+
																	2*fy31*gridded_density[j+1][k]+
																	2*fy12*gridded_density[j][k+1]+
																	fy12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fy31*gridded_density[j+1][k]+
																		2*fy31*gridded_density[j][k]+
																		2*fy12*gridded_density[j+1][k+1]+
																		fy12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fy13*gridded_density[j+1][k+1]+
																			2*fy12*gridded_density[j+1][k]+
																			2*fy13*gridded_density[j][k+1]+
																			fy12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fy13*gridded_density[j][k+1]+
																		2*fy13*gridded_density[j+1][k+1]+
																		2*fy12*gridded_density[j][k]+
																		fy12*gridded_density[j+1][k])/36*dV;	

					}
					else if (j==((int)p2Grid.e1()-1)){  //right edge, no corner
						dV = grid->cellVolume(j,k);
						fx = (int)p2Grid.e1()-p2Grid.e1()+1;

						fx13 = cube(1-fx);
						fx31 = 1-cube(fx);
						fx12 = (1+2*fx)*sqr(1-fx);

						int_density[j][k] += (4*fx31*gridded_density[j][k]+
																	2*fx12*gridded_density[j+1][k]+
																	2*fx31*gridded_density[j][k+1]+
																	fx12*gridded_density[j+1][k+1])/36*dV;
						int_density[j+1][k] += (4*fx13*gridded_density[j+1][k]+
																		2*fx12*gridded_density[j][k]+
																		2*fx13*gridded_density[j+1][k+1]+
																		fx12*gridded_density[j][k+1])/36*dV;
						int_density[j+1][k+1] += (4*fx13*gridded_density[j+1][k+1]+
																			2*fx13*gridded_density[j+1][k]+
																			2*fx12*gridded_density[j][k+1]+
																			fx12*gridded_density[j][k])/36*dV;
						int_density[j][k+1] += (4*fx31*gridded_density[j][k+1]+
																		2*fx12*gridded_density[j+1][k+1]+
																		2*fx31*gridded_density[j][k]+
																		fx12*gridded_density[j+1][k])/36*dV;	
					}
				}
				else{
					dV = grid->cellVolume(j,k);
					int_density[j][k] += (4*gridded_density[j][k]+2*gridded_density[j+1][k]+
																2*gridded_density[j][k+1]+gridded_density[j+1][k+1])/36*dV;
					int_density[j+1][k] += (4*gridded_density[j+1][k]+2*gridded_density[j][k]+
																	2*gridded_density[j+1][k+1]+gridded_density[j][k+1])/36*dV;
					int_density[j+1][k+1] += (4*gridded_density[j+1][k+1]+2*gridded_density[j+1][k]+
																		2*gridded_density[j][k+1]+gridded_density[j][k])/36*dV;
					int_density[j][k+1] += (4*gridded_density[j][k+1]+2*gridded_density[j+1][k+1]+
																	2*gridded_density[j][k]+gridded_density[j+1][k])/36*dV;
				}										
			}
		}
	
	for (j=(int)p1Grid.e1(); j<=MIN((int)p2Grid.e1(),J); j++)
		for (k=(int)p1Grid.e2(); k<=MIN((int)p2Grid.e2(),K); k++){
			total_density += int_density[j][k]; // this is really total charge
			int_density[j][k] /= grid->get_halfCellVolumes()[j][k];
			if (grid->axis()) {
				if (k == 0) int_density[j][k] *= 0.5/grid->get_radial0();
				else if (k == K) int_density[j][k] *= grid->get_radialK();
			}
			fields->setloaddensity(species->getID(),j,k,int_density[j][k]);
			if (np2c==0)
				fields->setbgrho(j,k,q*int_density[j][k]);
			if (gridded_density[j][k]>max_density) 
			  {
			    max_density=gridded_density[j][k];
			  }
		}

	if (np2c&&(numberParticles ==0))
		numberParticles = (int)(total_density/np2c+0.5);
	if (numberParticles){
	  if(max_density < 1E-16) { max_density = 1; }
		for (j=(int)p1Grid.e1(); j<=MIN((int)p2Grid.e1(),J); j++)
		  {
		    for (k=(int)p1Grid.e2(); k<=MIN((int)p2Grid.e2(),K); k++)
		      {
			gridded_density[j][k] /= max_density;
		      }
		  }
	}
	
	for (j=0; j<=J; j++)
	  delete [] int_density[j];
	delete [] int_density;
	int_density = 0;
	
}


Scalar Load::gridded_load(Vector2& x)
{
	Vector2 gridx = grid->getGridCoords(x);
	return grid->interpolateBilinear(gridded_density,gridx);
}

void Load::init_evaluator()
{
	int j,k;
	Vector2 x;
	Vector2 p1Grid = grid->getGridCoords(p1);
	Vector2 p2Grid = grid->getGridCoords(p2)+Vector2(1,1);
	for (j=(int)p1Grid.e1(); j<=MIN((int)p2Grid.e1(),J); j++)
	  {
	    for (k=(int)p1Grid.e2(); k<=MIN((int)p2Grid.e2(),K); k++)
	      { 
	      x = grid->getMKS(j,k);
	      adv_eval->add_variable("x1",x.e1());
	      adv_eval->add_variable("x2",x.e2());
	      gridded_density[j][k] = adv_eval->Evaluate(analyticF.c_str());
	      }
	  }
}

// Compute the coefficients used later in r-z geometry as needed

void Load::set_coefficients()
{
   if (rz) {
      rMin = p1.e2();
      rMax = p2.e2();
      rMinSqr = rMin * rMin;
      drSqr   = rMax * rMax - rMinSqr;
   }
}

void Load::init_default()
{
   volume = deltaP.e1();
   if(rz) {
      volume *= M_PI*drSqr;
   }
   else 
		volume *=deltaP.e2();

   if (np2c)
		numberParticles = (int)(volume*density/np2c+0.5);

	int j,k;
	Vector2 gridp1, gridp2;
	gridp1 = grid->getGridCoords(p1);
	gridp2 = grid->getGridCoords(p2);
	Vector2 upper,lower;
	for (j=((int)(gridp1.e1()+.5)); j<=((int)(gridp2.e1()+.5)); j++)
		for (k=((int)(gridp1.e2()+.5)); k<=((int)(gridp2.e2()+.5)); k++)
			gridded_density[j][k] = density;

}