#ifndef _BOUNDARY_H
#define _BOUNDARY_H



inline double nwall(double xshift, double nlength, double rexit, double rinlet, double x, double y, double z) //FHE Nozzle Wall Equation
{
    double rslope = ((rexit-rinlet)/nlength) ;
    double fnwall = rinlet + rslope* abs(x-xshift)  - sqrt(y * y + z * z);
    return (fnwall);
}


/* */


template<typename T>
struct init_boundary_neighbors_wall4BC //KGK
{
	const IBC_data& IBC_dat ;
	
	init_boundary_neighbors_wall4BC( const IBC_data& IBC_Da ) : IBC_dat(IBC_Da) { }
	void operator()(T& lattice) const
	{
		const xyz_type<int>& idx = lattice.idx;
		double nwallval;
		
		double xx = double( idx.x ) ;
		double yy = double( lattice.center_y - idx.y ) ;
		double zz = double( lattice.center_z - idx.z ) ;
		
		
		if (IBC_dat.nlength >= idx.x)
		{
			nwallval = 10000.0;
			nwallval = nwall(IBC_dat.xshift, IBC_dat.nlength, IBC_dat.rexit, IBC_dat.rinlet, xx, yy, zz);
			
			if (nwallval >= - sqrt(2.0) && nwallval <= 0.0)
			{
				lattice.onBoundaryTag = 1 ;										//JCR: set tag if @ boundary/shape/wall; set after .._periodic_BC
			}
		}
	}
} ;


template<typename T>
struct init_boundary_neighbors_wf
{
  void operator()(T& lattice) const
	{
		const xyz_type<int>& idx = lattice.idx;
		
		if (lattice.is_on_face(Xn))		//true if idx.x == 0
		{
			if (lattice.is_velinlet())
			{
				const int y_offset[5] = { 0,  0,  0,  1, -1};
				const int z_offset[5] = { 0,  1, -1,  0,  0};
				
				for (int i=0; i<5; i++)
					lattice.neighbors[i] =
					&(*(T::begin_iterator + T::offset(1, idx.y+y_offset[i], idx.z+z_offset[i])));
			} else {
				const int y_offset[5] = {  0,  0,  0, -2,  2};
				const int z_offset[5] = {  0,  2, -2,  0,  0};
				
				for (int i=0; i<5; i++)
					lattice.neighbors[i] =
					&(*(T::begin_iterator + T::offset(2, idx.y+y_offset[i], idx.z+z_offset[i])));	
			}
		}
		
		if (lattice.is_on_face(Xp))
		//	lattice.neighbors[0] =
			lattice.neighbors[ Xp ] =  // Shouldn't this be Xp if question is about Xp?
			&(*(T::begin_iterator + T::offset(idx.x-1, idx.y, idx.z)));
		
		const int y_offset[4] = {T::dims.y-2,        1,          idx.y,    idx.y};
		const int z_offset[4] = {      idx.z,    idx.z,    T::dims.z-2,        1};
		
		for (int i=Yn; i<=Zp; ++i)
		{
			if (lattice.is_on_face(i))
			//	lattice.neighbors[0] =
				lattice.neighbors[ i ] =   // Why did Nathan have 0 & not " i " or " i - 2 " here, as in _offset declared right before this loop & as used here, if looping over i? Strong fx on how strong RJ flow: lattice.neighbors[0]->fi used below
				&(*(T::begin_iterator + T::offset(idx.x, y_offset[i-2], z_offset[i-2])));
		}
	} 
};

template<typename T>
struct init_boundary_neighbors_periodic_BC //KGK
{
	const IBC_data& IBC_dat ;
	
	init_boundary_neighbors_periodic_BC( const IBC_data& IBC_Da ) : IBC_dat(IBC_Da) { }
	void operator()(T& lattice) const
  {
    const xyz_type<int>& idx = lattice.idx;
/*     
	  const int x_offset[6] = {T::dims.x-3,     2,        idx.x,  idx.x,       idx.x,  idx.x};
	  const int y_offset[6] = {idx.y,       idx.y,  T::dims.y-3,      2,       idx.y,  idx.y};   
	  const int z_offset[6] = {idx.z,       idx.z,        idx.z,  idx.z,  T::dims.z-3,     2};  
*/

      lattice.onBoundaryTag = 0;                                                // this must come 1st to ensure it is intialized to dflt

      if ( ( IBC_dat.periodicx*T::dims.x <= 3 ) || ( IBC_dat.periodicx*T::dims.y <= 3 ) || ( IBC_dat.periodicx*T::dims.z <= 3 ) ) return ;

	  const int x_offset[6] = {IBC_dat.periodicx*T::dims.x-3, IBC_dat.periodicx*2,      IBC_dat.periodicy*idx.x,	  IBC_dat.periodicy*idx.x,IBC_dat.periodicz*idx.x,		IBC_dat.periodicz*idx.x};
	  const int y_offset[6] = {IBC_dat.periodicx*idx.y,       IBC_dat.periodicx*idx.y,  IBC_dat.periodicy*T::dims.y-3,IBC_dat.periodicy*2,    IBC_dat.periodicz*idx.y,		IBC_dat.periodicz*idx.y};
	  const int z_offset[6] = {IBC_dat.periodicx*idx.z,       IBC_dat.periodicx*idx.z,  IBC_dat.periodicy*idx.z,	  IBC_dat.periodicy*idx.z,IBC_dat.periodicz*T::dims.z-3,IBC_dat.periodicz*2};

    for (int i=Xn; i<=Zp; ++i)  
    {
      if (lattice.is_on_face(i))
	  {
		lattice.neighbors[i] = &(*(T::begin_iterator + T::offset(x_offset[i], y_offset[i], z_offset[i])));
	    //KGK Store the indices for the neighbors for positive identification/checking
		lattice.neighborID[i][0] = lattice.neighbors[i]->mm; 
		lattice.neighborID[i][1] = lattice.neighbors[i]->nn; 
		lattice.neighborID[i][2] = lattice.neighbors[i]->rr;
		//KGK also set onBoundaryTag to 1
		lattice.onBoundaryTag = 1;
	  }
    }
  } 
};


template<typename T>
struct update_boundary_element_wf //JCR: 0 Order Accurate Extrapolation? E.g., last x face or right face using EXTRAPOLATION BOUNDARY CONDITIONS
{
  void operator()(T& lattice) const
  {
    std::copy(lattice.neighbors[0]->fi, lattice.neighbors[0]->fi + en , lattice.fi);
  }
};

template<typename T>
struct update_xn_boundary_element_wf //JCR: using "bounce-back" BC to pass a value of moving wall (via u @ specified location) to f
{
	const int IBc ;
	
	update_xn_boundary_element_wf(const int ibC) : IBc(ibC) { }
	
	void operator()(T& lattice) 
	{
		
#ifdef DEBUG
		for (int i=1; i<en; i++) if( (fabs(lattice.fi[i])>1.0e10) || isnan(lattice.fi[i]) || isinf(lattice.old_fi[i]) )
		{
			cout<<"update_xn_boundary_element_wf 1: fi[i]= "<<lattice.fi[i]<<" & old_fi="<<lattice.old_fi[i]<<" @ "<<lattice.idx<<" @ i="<<i<<endl;
            Pause( "update_xn_boundary_element_wf 1: boundary.h: 154: lattice.fi & old_fi" ) ;
		}
#endif
		
		if (lattice.is_velinlet()) 
			{
				double p[en];
				p[ 1] = cix[ 1] * lattice.u.x + ciy[ 1] * lattice.u.y + ciz[ 1] * lattice.u.z;
				p[11] = cix[11] * lattice.u.x + ciy[11] * lattice.u.y + ciz[11] * lattice.u.z;
				p[13] = cix[13] * lattice.u.x + ciy[13] * lattice.u.y + ciz[13] * lattice.u.z;
				p[ 7] = cix[ 7] * lattice.u.x + ciy[ 7] * lattice.u.y + ciz[ 7] * lattice.u.z;
				p[ 9] = cix[ 9] * lattice.u.x + ciy[ 9] * lattice.u.y + ciz[ 9] * lattice.u.z;
				
				lattice.fi[ 1] = lattice.neighbors[0]->fi[ 2] + 6.0 * lattice.rho * (1./18.0) * p[ 1];
				lattice.fi[11] = lattice.neighbors[1]->fi[14] + 6.0 * lattice.rho * (1./36.0) * p[11];
				lattice.fi[13] = lattice.neighbors[2]->fi[12] + 6.0 * lattice.rho * (1./36.0) * p[13];
				lattice.fi[ 7] = lattice.neighbors[3]->fi[10] + 6.0 * lattice.rho * (1./36.0) * p[ 7];
				lattice.fi[ 9] = lattice.neighbors[4]->fi[ 8] + 6.0 * lattice.rho * (1./36.0) * p[ 9];
			}
			else
			{
				//double sigma = 1.0 ; //lattice.T ; // but if Temp not updated yet (w/o a lattice.old_T or ... maybe later)
				//double Konst = lattice.rho * pow(2.0*pi*sigma,-1.5) ; //density-based 3D Maxwellian; L2_norm's lower w/f=const:steady-state
                if ( IBc == 7 )
				for ( int n = 0 ; n <= en - 1 ; n++ )
				{ 
					/*double	c_1 = cix[ n ] / ( pi ) , c_2 = ciy[ n ] / ( pi ) , c_3 = ciz[ n ] / ( pi ) ;
					 lattice.fi[n] = Fluid_Species::f_Eq(c_1,c_2,c_3,lattice.u.x,lattice.u.y,lattice.u.z,sigma,Konst) ; //3D Maxwellian: needs Temp
					 lattice.fi[n] = lattice.neighbors[ lattice.neighborID[n][0] ]->fi[n] ; //0th order extrapolation in x: [n][0=x] neighbor
					 lattice.fi[n] = lattice.neighbors[ lattice.neighborID[n][1] ]->fi[n] ; //0th order extrapolation in y: [n][1=y] neighbor
					 lattice.fi[n] = lattice.neighbors[ lattice.neighborID[n][2] ]->fi[n] ; //0th order extrapolation in z: [n][2=z] neighbor
					 lattice.fi[n] = 0;//Eliasson,B"Outflow boundary conditions for the Fourier transformed three-dimensional Vlasov-Maxwell system"JCP225(2007)1508-1532
					 lattice.fi[n] = lattice.neighbors[Xn]->fi[n] ; */ //bounce-back vs. specular in x? 
					lattice.fi[n] = lattice.old_fi[n] ; //f=const (steady-state flow OK?!); best L2_norm's
				} 
			} 
		
#ifdef DEBUG
		//std::cout<<"update_xn_boundary_element_wf 2: lattice.fi["<<0<<"]="<<lattice.fi[0]<<"\n";
		for ( int n = 0 ; n <= en - 1 ; n++ )
			if( isnan(lattice.fi[n]) || isinf(lattice.old_fi[n]) )
			{
				cout<<"update_xn_boundary_element_wf 2: lattice.neighborID[n][0]="<<lattice.neighborID[n][0]<<"\n" ;
				cout<<"update_xn_boundary_element_wf 3: fi="<<lattice.fi[n]<<" old_fi="<<lattice.old_fi[n]<<" T="<<lattice.T<<" @ "<<lattice.idx<<" @ n="<<n<<endl;
                Pause( "update_xn_boundary_element_wf: boundary.h: 198: lattice.neighborID & fi" ) ;
			}
#endif  
		
	}
};


template<typename T>
struct boundary_wall_x_f_wf //JCR: bounce-back off non-moving, non-flexing, non-permeable wall; NOT specular reflection
{
	const IBC_data& IBC_Dat ;
	
	boundary_wall_x_f_wf( const IBC_data& IBC_Da ) : IBC_Dat(IBC_Da) { }
	
	void operator()(T& lattice)
  {
    if (!lattice.is_velinlet())
		{
			lattice.fi[ 1] = lattice.neighbors[0]->fi[ 2]; 
			lattice.fi[11] = lattice.neighbors[1]->fi[14]; 
			lattice.fi[13] = lattice.neighbors[2]->fi[12]; 
			lattice.fi[ 9] = lattice.neighbors[3]->fi[ 8];
			lattice.fi[ 7] = lattice.neighbors[4]->fi[10];   
		}
		else
		{
			//double sigma = 1.0 ; //lattice.T ; // or a T_wall
			//double Konst = lattice.rho*pow(2.0*pi*sigma,-1.5);//http://research.nianet.org/~luo/Reprints-luo/Notes/TUBraunschweig_03/TUBraunschweig_03.pdf
            if ( IBC_Dat.IC_BC == 7 )
			for ( int n = 0 ; n <= en - 1 ; n++ )
			{ //3D Maxwellian: 
				//double	  c_1 = cix[ n ] / ( pi ) , c_2 = ciy[ n ] / ( pi ) , c_3 = ciz[ n ] / ( pi ) ;
				//lattice.fi[n] = Fluid_Species::f_Eq(c_1,c_2,c_3,0,0,0,sigma,Konst);//cout<<" fi["<<n<<"]="<<lattice.fi[n];
				lattice.fi[n] = lattice.old_fi[n] ;//if Temp not updated yet; L2_norm's lower w/f=const: steady-state
			}
		} 
	  
#ifdef DEBUG
	  for ( int n = 0 ; n <= en - 1 ; n++ )
	  if( isnan(lattice.fi[n]) || isinf(lattice.old_fi[n]) )
	  {
		  cout<<"boundary_wall_x_f_wf:fi="<<lattice.fi[n]<<" old_fi="<<lattice.old_fi[n]<<" T="<<lattice.T<<" @ "<<lattice.idx<<" n="<<n<<endl;
          Pause( "boundary_wall_x_f_wf: boundary.h: 242: lattice.old_fi, fi & T" ) ;
	  }
#endif  
	  
  }
};

/*
template<typename T>
struct nozzle_wall_bc //FHE: bounce-back off Nozzle Wall
{
    void operator()(T& lattice) const
    {

                                       //1  2   3  4   5  6   7   8   9 10   11  12  13 14  15  16  17 18
        const int x_offset[ en - 1 ] = {-1, 1,  0, 0,  0, 0, -1,  1, -1, 1,  -1,  1, -1, 1,  0,  0,  0, 0};
        const int y_offset[ en - 1 ] = { 0, 0, -1, 1,  0, 0, -1, -1,  1, 1,   0,  0,  0, 0, -1,  1, -1, 1};
        const int z_offset[ en - 1 ] = { 0, 0,  0, 0, -1, 1,  0,  0,  0, 0,  -1, -1,  1, 1, -1, -1,  1, 1};

        const xyz_type<int>& idx = lattice.idx;
        const xyz_type<int>& dims = lattice.dims;

        for (int i = 0; i < en - 1; i++)
        {
            //KGK: For boundary lattices, definition of neighbors depends on the boundary conditions
            //KGK: This is taken care of in the boundary functions. Hence, set dummies (self) for boundary lattices
            //KGK: gcc is forgiving (not crashing) but this change is required for visual C++
            if ((idx.x + x_offset[i]) >= 0 && (idx.y + y_offset[i]) >= 0 && (idx.z + z_offset[i]) >= 0 &&
                    (idx.x + x_offset[i]) < dims.x && (idx.y + y_offset[i]) < dims.y && (idx.z + z_offset[i]) < dims.z)

                lattice.neighbors[i] = &(*(T::begin_iterator + T::offset(idx.x + x_offset[i], idx.y + y_offset[i], idx.z + z_offset[i])));
            else
                lattice.neighbors[i] = &(*(T::begin_iterator));
            //KGK Store the indices for the neighbors for positive identification/checking
            lattice.neighborID[i][0] = lattice.neighbors[i]->mm;
            lattice.neighborID[i][1] = lattice.neighbors[i]->nn;
            lattice.neighborID[i][2] = lattice.neighbors[i]->rr;

        }


        xyz_dt& u          = lattice.u;
        const PDF_type& fi = lattice.fi;

        double nwallval;

        double xx = idx.x;
        double yy = T::center_y - (lattice.idx.y);
        double zz = T::center_z - (lattice.idx.z);


        if (nlength >= idx.x)
        {
            nwallval = 10000;
            nwallval = nwall(xshift, nlength, rexit, rinlet, xx, yy, zz);

            if (nwallval >= -sqrt(2) && nwallval <= 0)
                {
                    lattice.fi[ 1] = lattice.neighbors[0]->fi[ 2];
                    lattice.fi[ 2] = lattice.neighbors[1]->fi[ 1];
                    
                    lattice.fi[ 3] = lattice.neighbors[2]->fi[ 4];
                    lattice.fi[ 4] = lattice.neighbors[3]->fi[ 3];
                    
                    lattice.fi[ 5] = lattice.neighbors[4]->fi[ 6];
                    lattice.fi[ 6] = lattice.neighbors[5]->fi[ 5];
                                        
                    lattice.fi[ 7] = lattice.neighbors[6]->fi[10];
                    lattice.fi[ 10] = lattice.neighbors[9]->fi[ 7];

                    lattice.fi[ 9] = lattice.neighbors[8]->fi[ 8];
                    lattice.fi[ 8] = lattice.neighbors[7]->fi[ 9];

                    lattice.fi[11] = lattice.neighbors[10]->fi[14];
                    lattice.fi[14] = lattice.neighbors[13]->fi[11];

                    lattice.fi[12] = lattice.neighbors[11]->fi[13];
                    lattice.fi[13] = lattice.neighbors[12]->fi[12];

                    lattice.fi[15] = lattice.neighbors[14]->fi[18];
                    lattice.fi[18] = lattice.neighbors[17]->fi[15];

                    lattice.fi[16] = lattice.neighbors[15]->fi[17];
                    lattice.fi[17] = lattice.neighbors[16]->fi[16];
                   
                    
                }           
            
        }
    }
};

*/




//KGK Set value for one lattice cell
void dirichlet_bc_set_values (grid_type& grid, const char fieldType, const int mm, const int nn, const int rr, const double *bcvalue) //KGK
{		
	lattice_type* dirich_lattice;
	
	dirich_lattice = &(*(grid.begin() + lattice_type::offset(mm-1, nn-1, rr-1)));
	switch(fieldType)
	{
		case('A'): dirich_lattice->A        = *bcvalue;                                         break; //Magnetic potential
		case('B'): dirich_lattice->B        = *bcvalue;                                         break; //Magnetic field
		case('E'): dirich_lattice->E        = *bcvalue;                                         break; //Electric field
		case('u'): dirich_lattice->u        = *bcvalue;                                         break; //Velocity field: JCR; not f
		case('f'): for ( int i = 0 ; i < en ; i++ ) dirich_lattice->fi[ i ] = bcvalue[ 0 ] ;    break; //JCR: f
		case('P'): dirich_lattice->Phi      = bcvalue[ 0 ] ;                                    break; //Electric potential
		case('R'): dirich_lattice->rho      = bcvalue[ 0 ] ;                                    break; //Density
		case('J'): dirich_lattice->J        = *bcvalue;                                         break; //Current density
		case('c'): dirich_lattice->curlJ    = *bcvalue;                                         break; //Curl of current density
		case('g'): dirich_lattice->gradrho  = *bcvalue;                                         break; //Density gradient
	}
};


//KGK Set value for non-corner boundary cell
struct periodic_bc_set_value_non_corner_cell //KGK
{
	 const char fieldType;
	 const int locationTag;
	 const IBC_data IBC_dat;

	 periodic_bc_set_value_non_corner_cell (const char& fieldType1, const int locationTag1, const IBC_data IBC_Dat): 
	 fieldType(fieldType1), locationTag(locationTag1), IBC_dat(IBC_Dat) {}

     void operator()(lattice_type& lattice) const
     {
		if (!lattice.onBoundaryTag) return;							//Skip if not boundary cell
		if (locationTag==0 || locationTag == 1)
			if (lattice.is_at_corner()) return;						//Skip if boundary cell is at corner
		if (locationTag ==-1 && !lattice.is_at_corner()) return;	//Skip if boundary cell is not at corner

		lattice_type* pbc_neighbor;
/*
		 const int x_offset[6] = {1, -1,  0,  0,  0,  0};
		 const int y_offset[6] = {0,  0,	 1, -1,  0,  0};   
		 const int z_offset[6] = {0,  0,  0,  0,  1, -1}; 
*/
		 const int x_offset[6] = {IBC_dat.periodicx, -IBC_dat.periodicx,  0,  0,  0,  0};
		 const int y_offset[6] = {0,  0,	 IBC_dat.periodicy, -IBC_dat.periodicy,  0,  0};
		 const int z_offset[6] = {0,  0,  0,  0,  IBC_dat.periodicz, -IBC_dat.periodicz};
/*
     if ( ( fieldType == 'f' ) || ( fieldType == 'u' ) ) //JCR: f or u: same as below!
     {
         if ( lattice.idx.x == 0 ) // Xn or i=0 face
         {
             pbc_neighbor = &(*(lattice_type::begin_iterator + lattice_type::offset(lattice.dims.x-2, lattice.idx.y, lattice.idx.z)));
         }
         else if ( lattice.idx.x == lattice.dims.x-1 ) // Xp or i=Nx face
         {
             pbc_neighbor = &(*(lattice_type::begin_iterator + lattice_type::offset(1, lattice.idx.y, lattice.idx.z)));
         } else if ( lattice.idx.y == 0 ) // Yn or j=0 face
         {
             pbc_neighbor = &(*(lattice_type::begin_iterator + lattice_type::offset(lattice.idx.x, lattice.dims.y-2, lattice.idx.z)));
         }
         else if ( lattice.idx.y == lattice.dims.y-1 ) // Yp or i=Ny face
         {
             pbc_neighbor = &(*(lattice_type::begin_iterator + lattice_type::offset(lattice.idx.x, 1, lattice.idx.z)));
         } 
         else if ( lattice.idx.z == 0 ) // Zn or k=0 face
         {
             pbc_neighbor = &(*(lattice_type::begin_iterator + lattice_type::offset(lattice.idx.x, lattice.idx.y, lattice.dims.z-2)));
         } 
         else if ( lattice.idx.z == lattice.dims.z-1 ) // Zp or i=Nz face
         {
             pbc_neighbor = &(*(lattice_type::begin_iterator + lattice_type::offset(lattice.idx.x, lattice.idx.y, 1)));
         } 
         std::copy( pbc_neighbor->fi , pbc_neighbor->fi + en , lattice.fi ) ;
         lattice.u = pbc_neighbor->u;
     }
     else */
		for (int i=Xn; i<=Zp; ++i)
		  if (lattice.is_on_face(i)) 
		  {			  
			  int mm = lattice.neighborID[i][0]-1;  int nn = lattice.neighborID[i][1]-1;  int rr = lattice.neighborID[i][2]-1;
			  pbc_neighbor = &(*(lattice_type::begin_iterator + lattice_type::offset(mm+x_offset[i], nn+y_offset[i], rr+z_offset[i])));
			  switch(fieldType)
				{
					case('A'): lattice.A        = pbc_neighbor->A;                                  break; //Magnetic potential
					case('B'): lattice.B        = pbc_neighbor->B;                                  break; //Magnetic field
					case('E'): lattice.E        = pbc_neighbor->E;                                  break; //Electric field
					case('u'): lattice.u        = pbc_neighbor->u;                                  break; //Velocity field: JCR; not f
					case('f'): std::copy( pbc_neighbor->fi , pbc_neighbor->fi + en , lattice.fi ) ; break; //JCR: f
					case('P'): lattice.Phi      = pbc_neighbor->Phi;                                break; //Electric Potential
					case('J'): lattice.J        = pbc_neighbor->J;                                  break; //Current Density
					case('c'): lattice.curlJ    = pbc_neighbor->curlJ;                              break; //curl of Current Density
					case('g'): lattice.gradrho  = pbc_neighbor->gradrho;                            break; //Gradient of Density
					case('R'): lattice.rho      = pbc_neighbor->rho;                                break; //Density
				}
/* #ifdef DEBUG
              cout<<"\nmm="<<mm<<" nn="<<nn<<" rr="<<rr<<" idx="<<lattice.idx<<" on face i="<<i;
              Pause( "periodic_bc_set_value_non_corner_cell: boundary.h: 436: periodic BC neighbor & idx" ) ;
#endif */
		  }
  }
};

void periodic_bc_set_values(grid_type& grid, const char fieldType, const IBC_data IBC_Dat) //KGK
{
	int locationTag;
	//Set values for all non-corner boundary cells
	locationTag = 1;
	std::for_each(grid.begin(), grid.end(), periodic_bc_set_value_non_corner_cell(fieldType, locationTag, IBC_Dat)); 

	//Set values for Zn,Zp cells and corner boundary cells
	locationTag = 0;
	std::for_each(grid.begin(), grid.end(), periodic_bc_set_value_non_corner_cell(fieldType, locationTag, IBC_Dat));

	//Set values for Zn,Zp cells and corner boundary cells
	locationTag = -1;
	std::for_each(grid.begin(), grid.end(), periodic_bc_set_value_non_corner_cell(fieldType, locationTag, IBC_Dat));
  
}

void boundary_conditions_f(grid_type& grid, const IBC_data IBC_Dat )
{
	const int pad_map[6] = {1, 1, 0, 0, 0, 0};
	const int BIC					 = IBC_Dat.IC_BC ;
	const int when					 = IBC_Dat.when ;
	const int when_to_start			 = IBC_Dat.when2begin ;
	const int when_to_stop			 = IBC_Dat.when2end ;
	const xyz_type<int> PhysicalGrid = IBC_Dat.Physical_Grid ;
	const int where_to_start		 = IBC_Dat.where.x ;
	xyz_type<int> loc				 = IBC_Dat.where ;
	double A_BC[3] = { 0.0 , 0.0 , 0.0 } ; // A default field BC 
	
	switch(IBC_Dat.IC_BC) //JCR: initialize flow & E&M fields
	{
		case 0 : // JCR: read data from result_all
			break;
		default :
		case 1 : //JCR Lattice Boltzman Method MRT & RPJ
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, pad_map[Xn]),
						  face_iterator<lattice_type>(grid.end()), 
						  update_xn_boundary_element_wf<lattice_type>( BIC ));	//only performs for velocity inlet

			for (int i=Xp; i<=Zp; ++i)
			{
				std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
							  face_iterator<lattice_type>(grid.end()),
							  update_boundary_element_wf<lattice_type>());		//updates info only; JCR: extrapolation
			}
			
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, 1),
						  face_iterator<lattice_type>(grid.end()),  
						  boundary_wall_x_f_wf<lattice_type>( IBC_Dat ));		//bounce-back off wall surrounding inlet
			break;
		case 2 : //JCR: turbulent velocity field
		case 3 : //JCR: Taylor-Green vortex IC; periodic flow patterns
		case 9 : //JCR: turbulent velocity & magnetic fields
			periodic_bc_set_values(grid, 'A' , IBC_Dat);
			periodic_bc_set_values(grid, 'B' , IBC_Dat);
			periodic_bc_set_values(grid, 'E' , IBC_Dat);
			periodic_bc_set_values(grid, 'P' , IBC_Dat);
			periodic_bc_set_values(grid, 'u' , IBC_Dat); 
			periodic_bc_set_values(grid, 'f' , IBC_Dat); 
			periodic_bc_set_values(grid, 'J' , IBC_Dat);
			periodic_bc_set_values(grid, 'c' , IBC_Dat);
			periodic_bc_set_values(grid, 'g' , IBC_Dat);
			periodic_bc_set_values(grid, 'R' , IBC_Dat);
            break;
		case 4 : //JCR: Uniform field(s) IC; inlet upstream; periodic or extrapolation & extrapolation downstream for freely flowing
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, pad_map[Xn]),
						  face_iterator<lattice_type>(grid.end()), 
						  update_xn_boundary_element_wf<lattice_type>( BIC ));	//only performs for velocity inlet
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xp, PhysicalGrid+2, pad_map[Xp]),
						  face_iterator<lattice_type>(grid.end()),
						  update_boundary_element_wf<lattice_type>());		//updates info only; JCR: extrapolation
			if ( IBC_Dat.periodicy == 1 )
			{
				periodic_bc_set_values(grid, 'A' , IBC_Dat);
				periodic_bc_set_values(grid, 'B' , IBC_Dat);
				periodic_bc_set_values(grid, 'E' , IBC_Dat);
                periodic_bc_set_values(grid, 'P' , IBC_Dat);
				periodic_bc_set_values(grid, 'u' , IBC_Dat);
				periodic_bc_set_values(grid, 'f' , IBC_Dat);
                periodic_bc_set_values(grid, 'J' , IBC_Dat);
                periodic_bc_set_values(grid, 'c' , IBC_Dat);
                periodic_bc_set_values(grid, 'g' , IBC_Dat);
                periodic_bc_set_values(grid, 'R' , IBC_Dat);
			}
			else
			{
				for (int i=Yn; i<=Yp; ++i)
				{
					std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
								  face_iterator<lattice_type>(grid.end()),
								  update_boundary_element_wf<lattice_type>());	//updates info only; JCR: extrapolation
				}
			}
			if ( IBC_Dat.periodicz == 1 )
			{
				periodic_bc_set_values(grid, 'A' , IBC_Dat);
				periodic_bc_set_values(grid, 'B' , IBC_Dat);
				periodic_bc_set_values(grid, 'E' , IBC_Dat);
                periodic_bc_set_values(grid, 'P' , IBC_Dat);
				periodic_bc_set_values(grid, 'u' , IBC_Dat);
				periodic_bc_set_values(grid, 'f' , IBC_Dat);
                periodic_bc_set_values(grid, 'J' , IBC_Dat);
                periodic_bc_set_values(grid, 'c' , IBC_Dat);
                periodic_bc_set_values(grid, 'g' , IBC_Dat);
                periodic_bc_set_values(grid, 'R' , IBC_Dat);
			}
			else
			{
				for (int i=Zn; i<=Zp; ++i)
				{
					std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
								  face_iterator<lattice_type>(grid.end()),
								  update_boundary_element_wf<lattice_type>());	//updates info only; JCR: extrapolation
				}
			}
			break;
	    case 5 : //JCR: Plane shock @ xShock IC; 
			periodic_bc_set_values(grid, 'A' , IBC_Dat);
			periodic_bc_set_values(grid, 'B' , IBC_Dat);
			periodic_bc_set_values(grid, 'E' , IBC_Dat);
			periodic_bc_set_values(grid, 'P' , IBC_Dat);
			periodic_bc_set_values(grid, 'J' , IBC_Dat);
			periodic_bc_set_values(grid, 'c' , IBC_Dat);
			periodic_bc_set_values(grid, 'g' , IBC_Dat);
			periodic_bc_set_values(grid, 'R' , IBC_Dat);
			if ( when > when_to_start ) 
			{
				A_BC[ 0 ] = IBC_Dat.u0.x * when/when_to_stop * 0.1 ; 
				A_BC[ 1 ] = IBC_Dat.u0.y ; 
				A_BC[ 2 ] = IBC_Dat.u0.z ;
				loc.x	  = where_to_start ;//+ int(A_BC[ 0 ] * when/when_to_stop);//location of moving shock front
				std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, pad_map[Xn]),
							  face_iterator<lattice_type>(grid.end()), 
							  update_xn_boundary_element_wf<lattice_type>(BIC));//only performs for velocity inlet
				for( int j = 1 ; j < PhysicalGrid.y + 1 ; j++ )
					for( int k = 1 ; k < PhysicalGrid.z + 1 ; k++ )
						dirichlet_bc_set_values(grid, 'u', loc.x, j, k, A_BC);	//Apply varying u_inlet(=u_Shock) bc @ an x & associated yz plane
			}
			for (int i=Xn; i<=Xp; ++i)
			{
				std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
							  face_iterator<lattice_type>(grid.end()),  
							  boundary_wall_x_f_wf<lattice_type>( IBC_Dat ));	//bounce-back off walls
			}
			for (int i=Yn; i<=Zp; ++i)
			{
				std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
							  face_iterator<lattice_type>(grid.end()),
							  update_boundary_element_wf<lattice_type>());		//updates info only; JCR: extrapolation
			}
			break;
		case 6 : // JCR: f = 0
			A_BC[ 0 ] = 0.0 ; 
            dirichlet_bc_set_values(grid, 'f', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC ) ; // Set f = 0
			break;
		case 7 : // JCR: f = constant
            for (int i=Xn; i<=Zp; ++i)
            {
                std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
                              face_iterator<lattice_type>(grid.end()),
                              update_boundary_element_wf<lattice_type>());      //updates info only; JCR:  f = constant = 0th order extrapolation
            }
			break;
		case 8 : //JCR: Supersonic flow field(s) IC; periodic approximates infinite extent 3D but " is a mist for hyperbolic M>1 flow
			A_BC[ 0 ] = IBC_Dat.B0.x ; 
			A_BC[ 1 ] = IBC_Dat.B0.y ; 
			A_BC[ 2 ] = IBC_Dat.B0.z ;
			dirichlet_bc_set_values(grid, 'B', loc.x, loc.y, loc.z, A_BC);		//Apply B0 bc @ an x & associated yz 
			A_BC[ 0 ] = IBC_Dat.E0.x ; 
			A_BC[ 1 ] = IBC_Dat.E0.y ; 
			A_BC[ 2 ] = IBC_Dat.E0.z ;
			dirichlet_bc_set_values(grid, 'E', loc.x, loc.y, loc.z, A_BC);		//Apply E0 bc @ an x & associated yz 
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, pad_map[Xn]),
						  face_iterator<lattice_type>(grid.end()), 
						  update_xn_boundary_element_wf<lattice_type>( BIC ));	//only performs for velocity inlet
			for (int i=Xp; i<=Zp; ++i)
			{
				std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
							  face_iterator<lattice_type>(grid.end()),
							  update_boundary_element_wf<lattice_type>());		//updates info only; JCR: extrapolation
			}
			break;
		case 10 ://JCR: inlet upstream & extrapolation downstream for free-flowing; walls elsewhere for channel flow; Hartmann, Couette, Poiseuille
			A_BC[ 0 ] = 0.0 ;  // dflt "grounded" wall
			A_BC[ 1 ] = 0.0 ; 
			A_BC[ 2 ] = 0.0 ;
            dirichlet_bc_set_values(grid, 'A', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: magnetic potential
            dirichlet_bc_set_values(grid, 'B', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: magnetic induction
            dirichlet_bc_set_values(grid, 'E', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: Electric field
            dirichlet_bc_set_values(grid, 'P', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: Electric potential
            dirichlet_bc_set_values(grid, 'J', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: current density
            dirichlet_bc_set_values(grid, 'c', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: curls of "
            dirichlet_bc_set_values(grid, 'g', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: grad(rho)
            A_BC[ 0 ] = 1.0 ; //Temporary: would only work for incompressible flow
            dirichlet_bc_set_values(grid, 'R', PhysicalGrid.x+1, PhysicalGrid.y+1, PhysicalGrid.z+1, A_BC);// Apply the bc values: Density
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, pad_map[Xn]),
						  face_iterator<lattice_type>(grid.end()), 
						  update_xn_boundary_element_wf<lattice_type>( BIC ));	//only performs for velocity inlet
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xp, PhysicalGrid+2, pad_map[Xp]),
						  face_iterator<lattice_type>(grid.end()),
						  update_boundary_element_wf<lattice_type>());          //updates info only; JCR: extrapolation
			for (int i=Yn; i<=Yp; ++i)
			{
				std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
							  face_iterator<lattice_type>(grid.end()),
							  boundary_wall_x_f_wf<lattice_type>( IBC_Dat ));   //only when NOT @ velocity inlet; bounce-back downstream?x_max?
			}
			for (int i=Zn; i<=Zp; ++i)
			{
				std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
							  face_iterator<lattice_type>(grid.end()),
							  boundary_wall_x_f_wf<lattice_type>( IBC_Dat ));   //bounce-back side walls
			}
			break;
		case 11 : //JCR Lattice Boltzman Method MRT & RPJ w/|| not necessarily conducting plates; if conducting must account differently!
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, pad_map[Xn]),
						  face_iterator<lattice_type>(grid.end()), 
						  update_xn_boundary_element_wf<lattice_type>( BIC ));	//only performs for velocity inlet

			if ( IBC_Dat.periodicy == 1 )
			{
				periodic_bc_set_values(grid, 'A' , IBC_Dat);
				periodic_bc_set_values(grid, 'B' , IBC_Dat);
				periodic_bc_set_values(grid, 'E' , IBC_Dat);
				periodic_bc_set_values(grid, 'u' , IBC_Dat);
				periodic_bc_set_values(grid, 'f' , IBC_Dat);
                periodic_bc_set_values(grid, 'P' , IBC_Dat);
                periodic_bc_set_values(grid, 'J' , IBC_Dat);
                periodic_bc_set_values(grid, 'c' , IBC_Dat);
                periodic_bc_set_values(grid, 'g' , IBC_Dat);
                periodic_bc_set_values(grid, 'R' , IBC_Dat);
			}
			else
			{
				for (int i=Yn; i<=Yp; ++i)
				{
					std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
								  face_iterator<lattice_type>(grid.end()),
								  update_boundary_element_wf<lattice_type>());	//updates info only; JCR: extrapolation
				}
			}
			if ( IBC_Dat.periodicz == 1 )
			{
				periodic_bc_set_values(grid, 'A' , IBC_Dat);
				periodic_bc_set_values(grid, 'B' , IBC_Dat);
				periodic_bc_set_values(grid, 'E' , IBC_Dat);
				periodic_bc_set_values(grid, 'u' , IBC_Dat);
				periodic_bc_set_values(grid, 'f' , IBC_Dat);
                periodic_bc_set_values(grid, 'P' , IBC_Dat);
                periodic_bc_set_values(grid, 'J' , IBC_Dat);
                periodic_bc_set_values(grid, 'c' , IBC_Dat);
                periodic_bc_set_values(grid, 'g' , IBC_Dat);
                periodic_bc_set_values(grid, 'R' , IBC_Dat);
			}
			else
			{
				for (int i=Zn; i<=Zp; ++i)
				{
					std::for_each(face_iterator<lattice_type>(grid.begin(), i, PhysicalGrid+2, pad_map[i]),
								  face_iterator<lattice_type>(grid.end()),
								  update_boundary_element_wf<lattice_type>());	//updates info only; JCR: extrapolation
				}
			}

			std::for_each(face_iterator<lattice_type>(grid.begin(), Xn, PhysicalGrid+2, 1),
						  face_iterator<lattice_type>(grid.end()),  
						  boundary_wall_x_f_wf<lattice_type>( IBC_Dat ));		//bounce-back off not necessarily conducting wall surrounding inlet
			std::for_each(face_iterator<lattice_type>(grid.begin(), Xp, PhysicalGrid+2, 1),
						  face_iterator<lattice_type>(grid.end()),  
						  boundary_wall_x_f_wf<lattice_type>( IBC_Dat ));		//bounce-back off not necessarily conducting wall opposite inlet
			break;
	}
	
	//(solver<3)
	//    std::for_each(grid.begin(), grid.end(), nozzle_wall_bc<lattice_type>());
	
};

#endif