Namespace containing helper functions for Stefan Boltzmann radiation. More...

Functions
Vector< double >	Max_coord (2,-DBL_MAX)
	Max bin coords. More...

Vector< double >	Min_coord (2, DBL_MAX)
	Min bin coords. More...

Vector< double >	Dx (2, 0.0)
	Increments in bin. More...

void	bin_helper (const Vector< Vector< double > > &ray_vertex, const bool &populate_bin, Vector< std::pair< unsigned, unsigned > > &intersected_bin, FiniteElement *el_pt, std::ofstream &outfile)

void	doc_bins (std::ofstream &bin_file)

void	doc_sample_points (std::ofstream &outfile, const Vector< FiniteElement * > &sb_face_element_pt)

void	setup_stefan_boltzmann_visibility (const Vector< FiniteElement * > &sb_face_element_pt)

Variables
Vector< Vector< std::set< FiniteElement * > > >	Element_in_bin
	Bin to store pointers to finite elements in. More...

unsigned	Nx_bin =1000
	Number of bins in x direction. More...

unsigned	Ny_bin =1000
	Number of bins in y direction. More...

unsigned	Nsample =10

Detailed Description

Namespace containing helper functions for Stefan Boltzmann radiation.

/////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////

Function Documentation

◆ bin_helper()

void oomph::StefanBoltzmannHelper::bin_helper	(	const Vector< Vector< double > > &	ray_vertex,
		const bool &	populate_bin,
		Vector< std::pair< unsigned, unsigned > > &	intersected_bin,
		FiniteElement *	el_pt,
		std::ofstream &	outfile
	)

Helper function to setup or use bin that is used to locate intersections of rays with boundaries.

Args: – Two vectors, ray_vertex[0] and ray_vertex[1] which define the end points of a ray. – Boolean populate_bin: Element_in_bin by associating all bins that are intersected by the ray with the element pointed to by el_pt. If false: Return vector of bin indices that are intersected by the ray (retrived from previously set up bin structure) – Output stream outfile; if open, we write the intersected bins to that file but only if populate_bin=false. Only used for debugging...

  {
   
   // Actually plot
   bool plot_it=false;
   if (outfile.is_open())
    {
     if (populate_bin)
      {
       //Issue a warning
       OomphLibWarning(
        "Not outputting while bin is being populated...\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
     else
      {
       plot_it=true;
      }
    }
   
   // Step along ray in increasing y direction
   Vector<double> lower_ray_vertex=ray_vertex[0];
   Vector<double> upper_ray_vertex=ray_vertex[1];
   if (ray_vertex[1][1]<ray_vertex[0][1])
    {
     lower_ray_vertex=ray_vertex[1];
     upper_ray_vertex=ray_vertex[0];
    }
   
   // Output
   if (plot_it)
    {
     outfile
      << "ZONE T=\"ray\"\n"
      <<ray_vertex[0][0] << " " << ray_vertex[0][1] << "\n"
      <<ray_vertex[1][0] << " " << ray_vertex[1][1] << "\n";
    }
   
   // Add bin that contains the start point of the ray
   unsigned ix_start=unsigned((lower_ray_vertex[0]-Min_coord[0])/
                              (Max_coord[0]-Min_coord[0])*
                              double(Nx_bin));
   unsigned iy_start=unsigned((lower_ray_vertex[1]-Min_coord[1])/
                              (Max_coord[1]-Min_coord[1])*
                              double(Ny_bin));
   if (populate_bin)
    {
     Element_in_bin[ix_start][iy_start].insert(el_pt);
    }
   else
    {
     intersected_bin.push_back(std::make_pair(ix_start,iy_start));
    }
   
   // Vertical ray?
   if (lower_ray_vertex[0]==upper_ray_vertex[0])
    {
     // Add all the bins above the initial one until end point
     unsigned iy_end=unsigned((upper_ray_vertex[1]-Min_coord[1])/
                              (Max_coord[1]-Min_coord[1])*
                              double(Ny_bin));
     for (unsigned i=iy_start+1;i<=iy_end;i++)
      {
       if (populate_bin)
        {
         Element_in_bin[ix_start][i].insert(el_pt);
        }
       else
        {
         intersected_bin.push_back(std::make_pair(ix_start,i));
        }
      }
    }
   // Non-vertical ray
   else
    {
     // Get the (finite) slope
     double slope=
      (upper_ray_vertex[1]-lower_ray_vertex[1])/
      (upper_ray_vertex[0]-lower_ray_vertex[0]);
     
     // Current bin level
     unsigned iy=iy_start;
     
     // Upper and lower bounds of current bin level
     double y_this_bin_min=Min_coord[1]+double(iy  )*Dx[1];
     double y_this_bin_max=Min_coord[1]+double(iy+1)*Dx[1];
     
     // Keep going until we've moved through all the
     // bins beyond the y-level of the ray's end point
     int unit_offset=1;
     while (y_this_bin_min<upper_ray_vertex[1])
      {
       // Work out x-coordinate of intersection of ray with
       // other boundary of its current y-bin level
       double x_intersect=Max_coord[0];
       if (lower_ray_vertex[0]>upper_ray_vertex[0])
        {
         x_intersect=Min_coord[0];
        }
       if (slope!=0.0)
        {
         x_intersect=lower_ray_vertex[0]+
          (y_this_bin_max-lower_ray_vertex[1])/slope;
        }
       
       // Cut it off if it goes outside the bin structure
       if (x_intersect>Max_coord[0]) x_intersect=Max_coord[0] ;
       if (x_intersect<Min_coord[0]) x_intersect=Min_coord[0] ;
       
       // What's the x-bin index of that point
       unsigned ix_intersect=unsigned((x_intersect-Min_coord[0])/
                                      (Max_coord[0]-Min_coord[0])*
                                      double(Nx_bin));
  
       // limits for x bins:
       unsigned i_lo=ix_start;
       unsigned i_hi=ix_intersect;
  
       // Swap if going to the left and apply 
       // unit offset. Equal to one initially because the first
       // bin has already been filled; in subsequent rows of
       // bins, we add all of them.
       if (i_lo>i_hi)
        {
         i_lo=ix_intersect;
         i_hi=ix_start-unit_offset;
        }
       else
        {
         i_lo+=unit_offset;
        }
       unit_offset=0;      
  
       // ...but don't fall off the end... 
       if (i_hi==Nx_bin) i_hi-=1;
  
       // Now add all the bins at this y-level
       for (unsigned i=i_lo;i<=i_hi;i++)
        {
         // .. but don't go beyond the end of the ray
         bool add_it=true;
         if (slope>0.0)
          {
           double x_left_end_bin =Min_coord[0]+double(i  )*Dx[0];
           if (x_left_end_bin>upper_ray_vertex[0])
            {
             add_it=false;
            }
          }
         else if (slope<0.0)
          {
           double x_right_end_bin=Min_coord[0]+double(i+1)*Dx[0];
           if (x_right_end_bin<upper_ray_vertex[0])
            {
             add_it=false;
            }
          }
         else
          {
           double x_left_end_bin =Min_coord[0]+double(i  )*Dx[0];
           double x_right_end_bin=Min_coord[0]+double(i+1)*Dx[0];
           if (x_left_end_bin>std::max(upper_ray_vertex[0],lower_ray_vertex[0]))
            {
             add_it=false;
            }
           else if (x_right_end_bin<std::min(upper_ray_vertex[0],
                                             lower_ray_vertex[0]))
            {
             add_it=false;
            }
          }
  
         if (add_it)
          {
           if (populate_bin)
            {
             Element_in_bin[i][iy].insert(el_pt);
            }
           else
            {
             intersected_bin.push_back(std::make_pair(i,iy));
            }
          }
        }
         
       // Now update the counters: x bin level starts at the
       // point of intersection with the next level...
       ix_start=ix_intersect;
         
       // ...upper and lower boundary moves one level up
       y_this_bin_min+=Dx[1];
       y_this_bin_max+=Dx[1];
         
       // ...and bump up the level itself
       iy++;
      }
    }
     
   // Plot the intersected bins?
   if (plot_it)
    {
     unsigned nbin=intersected_bin.size();
     for (unsigned i=0;i<nbin;i++)
      {
       unsigned ix=intersected_bin[i].first;
       unsigned iy=intersected_bin[i].second;
         
       double x_lo=Min_coord[0]+double(ix)*Dx[0];
       double x_hi=x_lo+Dx[0];
       double y_lo=Min_coord[1]+double(iy)*Dx[1];
       double y_hi=y_lo+Dx[1];
         
       outfile << "ZONE\n"
               << x_lo << " " << y_lo << "\n"
               << x_hi << " " << y_lo << "\n"
               << x_hi << " " << y_hi << "\n"
               << x_lo << " " << y_hi << "\n"
               << x_lo << " " << y_lo << "\n";
      }
    }
  }

References Dx(), Element_in_bin, i, max, Max_coord(), min, Min_coord(), Nx_bin, Ny_bin, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and GlobalFct::plot_it().

Referenced by setup_stefan_boltzmann_visibility().

◆ doc_bins()

void oomph::StefanBoltzmannHelper::doc_bins ( std::ofstream & bin_file )

  {
   // Loop over bins
   for (unsigned ix=0;ix<Nx_bin;ix++)
    {
     for (unsigned iy=0;iy<Ny_bin;iy++)
      {
       // Anybody at home?
       if (Element_in_bin[ix][iy].size()!=0)
        {
         double x_lo=Min_coord[0]+double(ix)*Dx[0];
         double x_hi=x_lo+Dx[0];
         double y_lo=Min_coord[1]+double(iy)*Dx[1];
         double y_hi=y_lo+Dx[1];
           
         bin_file << "ZONE I=2, J=2\n"
                  << x_lo << " " << y_lo << "\n"
                  << x_hi << " " << y_lo << "\n"
                  << x_lo << " " << y_hi << "\n"
                  << x_hi << " " << y_hi << "\n";
  
        }
      }
    }
  }

References Dx(), Element_in_bin, Min_coord(), Nx_bin, Ny_bin, and size.

Referenced by StefanBoltzmannProblem< ELEMENT >::setup_sb_radiation(), and StefanBoltzmannProblem< ELEMENT >::StefanBoltzmannProblem().

◆ doc_sample_points()

void oomph::StefanBoltzmannHelper::doc_sample_points	(	std::ofstream &	outfile,
		const Vector< FiniteElement * > &	sb_face_element_pt
	)

  {    
   // Vector for  coordinates of sample point
   Vector<double> sample_point(2);
   Vector<double> s(1);
  
   // Loop over all face elements
   unsigned nel=sb_face_element_pt.size();
   for (unsigned e=0;e<nel;e++)
    {
     StefanBoltzmannRadiationBase* el_pt=
      dynamic_cast<StefanBoltzmannRadiationBase*>(
       sb_face_element_pt[e]);
 #ifdef PARANOID 
     if (el_pt==0)
      {
       std::stringstream error_message;
       error_message << "Failed to cast possible intersecting element "
                     << e
                     << " to  StefanBoltzmannRadiationBase";
       throw OomphLibError(
        error_message.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
 #endif
       
     // Loop over sampling points
     for (unsigned j_sample=0;j_sample<Nsample;j_sample++)
      {
       // Get vector of local coordinates of plot point j
       // as second vertex
       el_pt->get_s_plot(j_sample,Nsample,s);
         
       // Get coordinate of vertex
       el_pt->interpolated_x(s,sample_point);
         
       outfile << sample_point[0] << " " << sample_point[1] << "\n";
      }
    }
  }

References e(), oomph::FiniteElement::get_s_plot(), oomph::FiniteElement::interpolated_x(), Nsample, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and s.

Referenced by StefanBoltzmannProblem< ELEMENT >::setup_sb_radiation(), and StefanBoltzmannProblem< ELEMENT >::StefanBoltzmannProblem().

◆ Dx()

Vector<double> oomph::StefanBoltzmannHelper::Dx	(	2	,
		0.	0
	)

Increments in bin.

Referenced by bin_helper(), doc_bins(), setup_stefan_boltzmann_visibility(), and Eigen::umfpack_get_numeric().

◆ Max_coord()

Vector<double> oomph::StefanBoltzmannHelper::Max_coord	(	2	,
		-	DBL_MAX
	)

Max bin coords.

Referenced by bin_helper(), and setup_stefan_boltzmann_visibility().

◆ Min_coord()

Vector<double> oomph::StefanBoltzmannHelper::Min_coord	(	2	,
		DBL_MAX
	)

Min bin coords.

Referenced by bin_helper(), doc_bins(), and setup_stefan_boltzmann_visibility().

◆ setup_stefan_boltzmann_visibility()

void oomph::StefanBoltzmannHelper::setup_stefan_boltzmann_visibility ( const Vector< FiniteElement * > & sb_face_element_pt )

Setup mutual visibility for Stefan Boltzmann radiation for all face elements (derived from StefanBoltzmannRadiationBase) contained in vector.

  {
   // Output file for debugging
   const bool plot_it=false;
   std::ofstream some_file;
   char filename[100];
  
   // Loop over all face elements to wipe previous information
   unsigned nel=sb_face_element_pt.size();
   for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
    {
     StefanBoltzmannRadiationBase* illuminated_el_pt=
      dynamic_cast<StefanBoltzmannRadiationBase*>(
       sb_face_element_pt[e_illuminated]);
 #ifdef PARANOID
     if (illuminated_el_pt==0)
      {
       std::stringstream error_message;
       error_message << "Failed to cast illuminated element "
                     << e_illuminated
                     << " to  StefanBoltzmannRadiationBase";
       throw OomphLibError(
        error_message.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
 #endif
      
     // Forget any previous information
     illuminated_el_pt->wipe_stefan_boltzmann_illumination_info();
    }
  
  
   double t_start=TimingHelpers::timer();
  
   // Vector to illuminated/ing Gauss points
   Vector<double> r_illuminated(2);
   Vector<double> s_illuminated(1);
   Vector<double> unit_normal_illuminated(2);
   Vector<double> r_illuminating(2);
   Vector<double> s_illuminating(1);
   Vector<double> unit_normal_illuminating(2);
  
   // Setup bin structure
   const bool use_bins=true;
   if (use_bins)
    {
     double t_start_bin=TimingHelpers::timer();
  
     // Each bin stores FiniteElements that have any sample point in it
     Element_in_bin.clear();
     Element_in_bin.resize(Nx_bin);
     for (unsigned i=0;i<Nx_bin;i++)
      {
       Element_in_bin[i].resize(Ny_bin);
      }
      
  
     // Reset max/min coords
     Max_coord.clear();
     Max_coord.resize(2,-DBL_MAX);
     Min_coord.clear();
     Min_coord.resize(2, DBL_MAX);
  
     // Initial sweep: Get max/min coords
     //----------------------------------
  
     // Loop over all face elements
     unsigned nel=sb_face_element_pt.size();
     for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
      {
       StefanBoltzmannRadiationBase* illuminated_el_pt=
        dynamic_cast<StefanBoltzmannRadiationBase*>(
         sb_face_element_pt[e_illuminated]);
 #ifdef PARANOID
       if (illuminated_el_pt==0)
        {
         std::stringstream error_message;
         error_message << "Failed to cast illuminated element "
                       << e_illuminated
                       << " to  StefanBoltzmannRadiationBase";
         throw OomphLibError(
          error_message.str(),
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
        }
 #endif
        
       // Loop over iluminated integration points
       unsigned nint=illuminated_el_pt->integral_pt()->nweight();
       for (unsigned ipt_illuminated=0;ipt_illuminated<nint;ipt_illuminated++)
        {
         // Local coordinate of integration point
         s_illuminated[0]=
          illuminated_el_pt->integral_pt()->knot(ipt_illuminated,0);
          
         // No recycling of shape fcts -- may as well call interpolated_x
         // directly
         illuminated_el_pt->interpolated_x(s_illuminated,r_illuminated);
          
         for (unsigned i=0;i<2;i++)
          {
           if (r_illuminated[i]>Max_coord[i]) Max_coord[i]=r_illuminated[i];
           if (r_illuminated[i]<Min_coord[i]) Min_coord[i]=r_illuminated[i];
          }
        }
      }
  
     // Allow for offset
     double percentage_offset=5.0;
     for (unsigned i=0;i<2;i++)
      {
       double dx=Max_coord[i]-Min_coord[i];
       Max_coord[i]=Max_coord[i]+percentage_offset/100.0*dx;
       Min_coord[i]=Min_coord[i]-percentage_offset/100.0*dx;
      }
  
     // Update increments
     Dx[0]=(Max_coord[0]-Min_coord[0])/double(Nx_bin);
     Dx[1]=(Max_coord[1]-Min_coord[1])/double(Ny_bin);
      
  
     const bool test_horizontal_and_vertical_rays=false;
     if (test_horizontal_and_vertical_rays)
      {
       
       // Test for left to right horizontal ray
       {
        // Populate bin by associating all bins intersected by
        // ray from one to the other sample vertex with the
        // element (note that the ray can span multiple bins!)
        Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
        bool populate_bin=true;
        std::ofstream dummy_file;
        Vector<Vector<double> > sample_vertex(2);
        sample_vertex[0].resize(2);
        sample_vertex[1].resize(2);
        sample_vertex[0][0]=-0.9;
        sample_vertex[0][1]=0.9;
        sample_vertex[1][0]=0.9;
        sample_vertex[1][1]=0.9;
        
        StefanBoltzmannRadiationBase* el_pt=
         dynamic_cast<StefanBoltzmannRadiationBase*>(
          sb_face_element_pt[0]);
        
        bin_helper(sample_vertex,
                   populate_bin,
                   dummy_intersected_bin,
                   el_pt,
                   dummy_file);
       }
       
       
       
       
       // Test for right to left horizontal ray
       {
        // Populate bin by associating all bins intersected by
        // ray from one to the other sample vertex with the
        // element (note that the ray can span multiple bins!)
        Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
        bool populate_bin=true;
        std::ofstream dummy_file;
        Vector<Vector<double> > sample_vertex(2);
        sample_vertex[0].resize(2);
        sample_vertex[1].resize(2);
        sample_vertex[0][0]= 0.9;
        sample_vertex[0][1]=-0.9;
        sample_vertex[1][0]=-0.9;
        sample_vertex[1][1]=-0.9;
        
        StefanBoltzmannRadiationBase* el_pt=
         dynamic_cast<StefanBoltzmannRadiationBase*>(
          sb_face_element_pt[0]);
        
        bin_helper(sample_vertex,
                   populate_bin,
                   dummy_intersected_bin,
                   el_pt,
                   dummy_file);
       }
       
       
       // Test for bottom to top vertical ray
       {
        // Populate bin by associating all bins intersected by
        // ray from one to the other sample vertex with the
        // element (note that the ray can span multiple bins!)
        Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
        bool populate_bin=true;
        std::ofstream dummy_file;
        Vector<Vector<double> > sample_vertex(2);
        sample_vertex[0].resize(2);
        sample_vertex[1].resize(2);
        sample_vertex[0][0]=-0.7;
        sample_vertex[0][1]=-0.9;
        sample_vertex[1][0]=-0.7;
        sample_vertex[1][1]= 0.9;
        
        StefanBoltzmannRadiationBase* el_pt=
         dynamic_cast<StefanBoltzmannRadiationBase*>(
          sb_face_element_pt[0]);
        
        bin_helper(sample_vertex,
                   populate_bin,
                   dummy_intersected_bin,
                   el_pt,
                   dummy_file);
       }
       
       
       // Test for to to bottom vertical ray
       {
        // Populate bin by associating all bins intersected by
        // ray from one to the other sample vertex with the
        // element (note that the ray can span multiple bins!)
        Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
        bool populate_bin=true;
        std::ofstream dummy_file;
        Vector<Vector<double> > sample_vertex(2);
        sample_vertex[0].resize(2);
        sample_vertex[1].resize(2);
        sample_vertex[0][0]= 0.7;
        sample_vertex[0][1]= 0.9;
        sample_vertex[1][0]= 0.7;
        sample_vertex[1][1]=-0.9;
        
        StefanBoltzmannRadiationBase* el_pt=
         dynamic_cast<StefanBoltzmannRadiationBase*>(
          sb_face_element_pt[0]);
        
        bin_helper(sample_vertex,
                   populate_bin,
                   dummy_intersected_bin,
                   el_pt,
                   dummy_file);
       }
       
      } // end of test for horizontal and vertical rays
  
     // Setup bin for checking intersections with rays
     //-----------------------------------------------
      
     // Loop over all face elements
     nel=sb_face_element_pt.size();
     for (unsigned e=0;e<nel;e++)
      {
       StefanBoltzmannRadiationBase* el_pt=
        dynamic_cast<StefanBoltzmannRadiationBase*>(
         sb_face_element_pt[e]);
 #ifdef PARANOID
       if (el_pt==0)
        {
         std::stringstream error_message;
         error_message << "Failed to cast possible intersecting element "
                       << e
                       << " to  StefanBoltzmannRadiationBase";
         throw OomphLibError(
          error_message.str(),
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
        }
 #endif
  
       // Loop over ALL sampling points along element
       Vector<Vector<double> > sample_vertex(2);
       sample_vertex[0].resize(2);
       sample_vertex[1].resize(2);
        
       // Get vector of local coordinates of plot point j
       // as first vertex
       Vector<double> s(1);
       unsigned j_sample=0;
       el_pt->get_s_plot(j_sample,Nsample,s);
        
       // Get coordinate of first vertex
       el_pt->interpolated_x(s,sample_vertex[0]);
        
       // Loop over remaining sampling points
       for (unsigned j_sample=1;j_sample<Nsample;j_sample++)
        {
         // Get vector of local coordinates of plot point j
         // as second vertex
         el_pt->get_s_plot(j_sample,Nsample,s);
          
         // Get coordinate of vertex
         el_pt->interpolated_x(s,sample_vertex[1]);
         
         // Populate bin by associating all bins intersected by
         // ray from one to the other sample vertex with the
         // element (note that the ray can span multiple bins!)
         Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
         bool populate_bin=true;
         std::ofstream dummy_file;
         bin_helper(sample_vertex,
                    populate_bin,
                    dummy_intersected_bin,
                    el_pt,
                    dummy_file);
          
         // Shift along
         sample_vertex[0]=sample_vertex[1];
        }
      }
      
     double t_end_bin=TimingHelpers::timer();
     oomph_info << "Time for setting up bin: "
                << t_end_bin-t_start_bin << std::endl;
    
     // Number of elements
     nel=sb_face_element_pt.size();
  
     // Assume all elements have the same number of integration points
     unsigned nintpt_all=sb_face_element_pt[0]->integral_pt()->nweight();
  
     // Exploit symmetry during setup
     const bool exploit_symmetry=true;
  
     // Helper lookup scheme to exploit symmetry of interaction:
     // If integration point i_ed in element e_ed is illuminatED by
     // integration point i_ing in element e_ing (which is the
     // illuminatING one!) then the reverse is true too
     Vector<Vector<Vector<Vector<unsigned> > > > aux;
     if (exploit_symmetry)
      {
       nintpt_all=sb_face_element_pt[0]->integral_pt()->nweight();
       aux.resize(nel);
       for (unsigned e=0;e<nel;e++)
        {
         aux[e].resize(nintpt_all);
         for (unsigned ipt=0;ipt<nintpt_all;ipt++)
          {
           aux[e][ipt].resize(nel);
          }
        }
      }
   
     // Loop over all face elements -- viewed as illuminated ones
     for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
      {
       StefanBoltzmannRadiationBase* illuminated_el_pt=
        dynamic_cast<StefanBoltzmannRadiationBase*>(
         sb_face_element_pt[e_illuminated]);
  
 #ifdef PARANOID
       if (illuminated_el_pt==0)
        {
         std::stringstream error_message;
         error_message << "Failed to cast illuminated element "
                       << e_illuminated
                       << " to  StefanBoltzmannRadiationBase";
         throw OomphLibError(
          error_message.str(),
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
        }
 #endif
      
       // Recast to FaceElement
       FaceElement* illuminated_face_el_pt=
        dynamic_cast<FaceElement*>(illuminated_el_pt);
      
       // Loop over iluminated integration points
       unsigned nint=illuminated_el_pt->integral_pt()->nweight();
       for (unsigned ipt_illuminated=0;ipt_illuminated<nint;ipt_illuminated++)
        {
         // Local coordinate of integration point
         s_illuminated[0]=
          illuminated_el_pt->integral_pt()->knot(ipt_illuminated,0);
        
         // No recycling of shape fcts -- may as well call interpolated_x
         // directly
         illuminated_el_pt->interpolated_x(s_illuminated,r_illuminated);
        
         // Get outer unit normal
         illuminated_face_el_pt->outer_unit_normal(s_illuminated,
                                                   unit_normal_illuminated);
        
         // Now loop over all other elements (the iluminating ones)
         // Note: this includes the current one because its integration
         // points may illuminate each other!
         unsigned e_lo=0;
         if (exploit_symmetry) e_lo=e_illuminated;
         for (unsigned e_illuminating=e_lo;e_illuminating<nel;e_illuminating++)
          {
           StefanBoltzmannRadiationBase* illuminating_el_pt=
            dynamic_cast<StefanBoltzmannRadiationBase*>(
             sb_face_element_pt[e_illuminating]);
  
 #ifdef PARANOID
           if (illuminating_el_pt==0)
            {
             std::stringstream error_message;
             error_message << "Failed to cast illuminating element "
                           << e_illuminating
                           << " to  StefanBoltzmannRadiationBase";
             throw OomphLibError(
              error_message.str(),
              OOMPH_CURRENT_FUNCTION,
              OOMPH_EXCEPTION_LOCATION);
            }
 #endif
          
           // Recast to FaceElement
           FaceElement* illuminating_face_el_pt=
            dynamic_cast<FaceElement*>(illuminating_el_pt);
          
           // Storage to accumulate indices of integration points in
           // illuminating face element that is visible from
           // current integration point
           Vector<unsigned> illuminating_integration_point_index;
          
           // Loop over iluminating integration points
           unsigned nint=illuminating_el_pt->integral_pt()->nweight();
           for (unsigned ipt_illuminating=0;ipt_illuminating<nint;
                ipt_illuminating++)
            {
             // Local coordinate of integration point
             s_illuminating[0]=
              illuminating_el_pt->integral_pt()->knot(ipt_illuminating,0);
            
             // No recycling of shape fcts -- may as well call interpolated_x
             // directly
             illuminating_el_pt->interpolated_x(s_illuminating,r_illuminating);
            
             // Get outer unit normal
             illuminating_face_el_pt->
              outer_unit_normal(s_illuminating,
                                unit_normal_illuminating);
             
             // Get vector from illuminating point to illuminated point
             Vector<double> ray(2);
             ray[0]=r_illuminated[0]-r_illuminating[0];
             ray[1]=r_illuminated[1]-r_illuminating[1];
            
  
             Vector<Vector<double> > ray_vertex(2);
             ray_vertex[0].resize(2);
             ray_vertex[1].resize(2);
             ray_vertex[0][0]=r_illuminated[0];
             ray_vertex[0][1]=r_illuminated[1];
             ray_vertex[1][0]=r_illuminating[0];
             ray_vertex[1][1]=r_illuminating[1];
  
             // Do we radiate away from the positive face of the
             // illuminating point?
             double dot_illuminating=
              (unit_normal_illuminating[0]*ray[0]+
               unit_normal_illuminating[1]*ray[1]);
             if (dot_illuminating>0.0)
              {
              
               // Do we radiate onto from the positive face of the
               // illuminated point?
               double dot_illuminated=
                (unit_normal_illuminated[0]*ray[0]+
                 unit_normal_illuminated[1]*ray[1]);
               if (dot_illuminated<0.0)
                {
  
                 // Does the ray (a finite-length segment) from
                 // radiating to radiated point intersect any other elements?
                               
                 // Vector of pairs of bin coordinates that
                 // intersect ray
                 Vector<std::pair<unsigned,unsigned> > intersected_bin;
                 if (use_bins)
                  {
                   if (plot_it)
                    {
                     sprintf(filename,"RESLT/latest_ray.dat");
                     some_file.open(filename);
                    }
  
                   // Find bins that are intersected by ray
                   Vector<std::pair<unsigned,unsigned> > intersected_bin;
                   bool populate_bin=false;
                   FiniteElement* dummy_el_pt=0;
                   bin_helper(ray_vertex,
                              populate_bin,
                              intersected_bin,
                              dummy_el_pt,
                              some_file);
                  
                   // End plot
                   if (plot_it)
                    {
                     some_file.close();
                     pause("done latest ray");
                    }
  
                   // Storage for vertices of possible intersection with ray
                   Vector<Vector<double> > segment_vertex(2);
                   segment_vertex[0].resize(2);
                   segment_vertex[1].resize(2);
  
                   // Search through all elements in bins along ray
                   bool have_intersection=false;
                   unsigned nbin=intersected_bin.size();
                   for (unsigned b=0;b<nbin;b++)
                    {
                     // Get bin indices
                     unsigned i=intersected_bin[b].first;
                     unsigned j=intersected_bin[b].second;
  
                     // Loop over elements in that bin
                     for (std::set<FiniteElement*>::iterator it=
                           Element_in_bin[i][j].begin();
                          it!=Element_in_bin[i][j].end();it++)
                      {
                       // Get possibly blocking element
                       StefanBoltzmannRadiationBase* el_pt=
                        dynamic_cast<StefanBoltzmannRadiationBase*>(*it);
                      
 #ifdef PARANOID
                       if (illuminated_el_pt==0)
                        {
                         std::stringstream error_message;
                         error_message
                          << "Failed to cast possibly blocking element "
                          << " to  StefanBoltzmannRadiationBase";
                         throw OomphLibError(
                          error_message.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
                        }
 #endif
                       // Skip intersections with illuming/illuminated element
                       if ((el_pt!=illuminated_el_pt)&&
                           (el_pt!=illuminating_el_pt))
                        {
                         // Loop over sampling points along element, only
                         // up to Nsample-1 because we're forming segment
                         // with current and next sampling point.
                         for (unsigned j_sample=0;j_sample<Nsample-1;j_sample++)
                          {
                           // Get cector of local coordinates of plot point j
                           // as first vertex
                           Vector<double> s(1);
                           el_pt->get_s_plot(j_sample,Nsample,s);
                          
                           // Get coordinate of first vertex
                           el_pt->interpolated_x(s,segment_vertex[0]);
                          
                           // Get cector of local coordinates of plot point j+1
                           // as second vertex
                           el_pt->get_s_plot(j_sample+1,Nsample,s);
                          
                           // Get coordinate of second vertex
                           el_pt->interpolated_x(s,segment_vertex[1]);
                          
                           // Check intersection
                           bool intersection=IntersectionChecker::intersects(
                            segment_vertex,ray_vertex);
                          
                           // Bail out
                           if (intersection)
                            {
                             have_intersection=true;
                             break;
                            }
                          } // end of loop over sampling points
                        } // endif for self-intersection
                       if (have_intersection) break;
                      } // end of loop over elements in bin
                     if (have_intersection) break;
                    }// End of loop over bins that contain ray
                  
                  
                   // No intersection
                   if (!have_intersection)
                    {
                     // Visible, so add integration point
                     illuminating_integration_point_index.
                      push_back(ipt_illuminating);
                    }
                  }
                 // No bins: Brute force search loop
                 else
                  {
                   // Intersection for star-shaped (non-convex) polygon
                   // can only be detected in O(n^2) operations.
                   Vector<Vector<double> > segment_vertex(2);
                   segment_vertex[0].resize(2);
                   segment_vertex[1].resize(2);
                  
                   // Loop over all segments to test for intersection
                   bool have_intersection=false;
                   for (unsigned e_intersect=0;e_intersect<nel;
                        e_intersect++)
                    {
                    
                     // Skip intersection with illuminating/ed elements
                     if (! ( (e_intersect==e_illuminated) ||
                             (e_intersect==e_illuminating) ) )
                      {
                      
                       // Loop over sampling points along element
                       for (unsigned j=0;j<Nsample-1;j++)
                        {
                         // Get cector of local coordinates of plot point j
                         // as first vertex
                         Vector<double> s(1);
                         sb_face_element_pt[e_intersect]->
                          get_s_plot(j,Nsample,s);
                        
                         // Get coordinate of first vertex
                         sb_face_element_pt[e_intersect]->
                          interpolated_x(s,segment_vertex[0]);
                        
                         // Get cector of local coordinates of plot point j+1
                         // as second vertex
                         sb_face_element_pt[e_intersect]->
                          get_s_plot(j+1,Nsample,s);
                        
                         // Get coordinate of second vertex
                         sb_face_element_pt[e_intersect]->
                          interpolated_x(s,segment_vertex[1]);
                        
                         // Check intersection
                         bool intersection=IntersectionChecker::intersects(
                          segment_vertex,ray_vertex);
                        
                         if (intersection)
                          {
                           have_intersection=true;
                           break;
                          }
                        }
                      }
                     if (have_intersection)
                      {
                       break;
                      }
                    }
                  
                   // No intersection
                   if (!have_intersection)
                    {
                     // Visible, so add integration point
                     illuminating_integration_point_index.
                      push_back(ipt_illuminating);
                    }
  
                  } // end if for brute force (rather than bin-based) search loop
  
                }
              }
            }
          
           // Done all possibly illuminating integration points in
           // current possibly illuminating element
           unsigned npt=illuminating_integration_point_index.size();
           if (npt>0)
            {
             // Add info
             illuminated_el_pt->add_stefan_boltzmann_illumination_info(
              ipt_illuminated,illuminating_el_pt,
              illuminating_integration_point_index);
  
             // Set up reverse scheme
             if (exploit_symmetry)
              {
               for (unsigned i_ing=0;i_ing<npt;i_ing++)
                {
                 aux[e_illuminating]
                  [illuminating_integration_point_index[i_ing]]
                  [e_illuminated].push_back(ipt_illuminated);
                }
              }
            }
          }
        }
      }
  
  
     // Now do other half of dependencies
     //----------------------------------
     if (exploit_symmetry)
      {
       // Loop over all face elements -- viewed as illuminated ones
       for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
        {
         StefanBoltzmannRadiationBase* illuminated_el_pt=
          dynamic_cast<StefanBoltzmannRadiationBase*>(
           sb_face_element_pt[e_illuminated]);
  
 #ifdef PARANOID
         if (illuminated_el_pt==0)
          {
           std::stringstream error_message;
           error_message << "Failed to cast illuminated element "
                         << e_illuminated
                         << " to  StefanBoltzmannRadiationBase";
           throw OomphLibError(
            error_message.str(),
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
          }
 #endif
  
         // Loop over iluminated integration points
         unsigned nint=illuminated_el_pt->integral_pt()->nweight();
         for (unsigned ipt_illuminated=0;ipt_illuminated<nint;ipt_illuminated++)
          {
           // Now loop over all other elements (the iluminating ones)
           for (unsigned e_illuminating=0;
                e_illuminating<e_illuminated;e_illuminating++)
            {
             StefanBoltzmannRadiationBase* illuminating_el_pt=
              dynamic_cast<StefanBoltzmannRadiationBase*>(
               sb_face_element_pt[e_illuminating]);
            
 #ifdef PARANOID
             if (illuminating_el_pt==0)
              {
               std::stringstream error_message;
               error_message << "Failed to cast illuminating element "
                             << e_illuminating
                             << " to  StefanBoltzmannRadiationBase";
               throw OomphLibError(
                error_message.str(),
                OOMPH_CURRENT_FUNCTION,
                OOMPH_EXCEPTION_LOCATION);
              }
 #endif
             // Get illuminating integration points
             Vector<unsigned> illuminating_integration_point_index=
              aux[e_illuminated][ipt_illuminated][e_illuminating];
             unsigned npt=illuminating_integration_point_index.size();
             if (npt>0)
              {
               // Add info
               illuminated_el_pt->add_stefan_boltzmann_illumination_info(
                ipt_illuminated,illuminating_el_pt,
                illuminating_integration_point_index);
              }
            }
          }
        }
      }
  
     double t_end=TimingHelpers::timer();
     oomph_info << "Time for setting up mutual Stefan Boltzmann radiation: "
                << t_end-t_start << std::endl;
    }
  
  }