temporary_stefan_boltzmann_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 // Header file for elements that are used to impose Stefan Boltzmann radiation
27 // on the surface of unsteady heat elements
28 
29 #ifndef OOMPH_hierher_HEAT_TRANSFER_AND_MELT_ELEMENTS_HEADER
30 #define OOMPH_hierher_HEAT_TRANSFER_AND_MELT_ELEMENTS_HEADER
31 
32 // Config header generated by autoconfig
33 #ifdef HAVE_CONFIG_H
34  #include <oomph-lib-config.h>
35 #endif
36 
37 //Standard libray headers
38 #include <cmath>
39 
40 // oomph-lib includes
41 #include "generic.h" // ../generic/Qelements.h"
42 
43 
44 namespace oomph
45 {
46 
47 //=======start_namespace==========================================
49 //================================================================
50 namespace IntersectionChecker
51 {
52 
57  bool intersects(const Vector<Vector<double> >& first_segment_vertex,
58  const Vector<Vector<double> >& second_segment_vertex,
59  const double& epsilon_parallel=1.0e-15)
60  {
61 
62  double denom = ((first_segment_vertex[1][1] - first_segment_vertex[0][1])*
63  (second_segment_vertex[1][0] - second_segment_vertex[0][0])) -
64  ((first_segment_vertex[1][0] - first_segment_vertex[0][0])*
65  (second_segment_vertex[1][1] - second_segment_vertex[0][1]));
66 
67  double nume_a = ((first_segment_vertex[1][0] - first_segment_vertex[0][0])*
68  (second_segment_vertex[0][1] - first_segment_vertex[0][1])) -
69  ((first_segment_vertex[1][1] - first_segment_vertex[0][1])*
70  (second_segment_vertex[0][0] - first_segment_vertex[0][0]));
71 
72  double nume_b = ((second_segment_vertex[1][0] - second_segment_vertex[0][0])*
73  (second_segment_vertex[0][1] - first_segment_vertex[0][1])) -
74  ((second_segment_vertex[1][1] - second_segment_vertex[0][1])*
75  (second_segment_vertex[0][0] - first_segment_vertex[0][0]));
76 
77  if(std::fabs(denom) < epsilon_parallel)
78  {
79  if( (std::fabs(nume_a) < epsilon_parallel) &&
80  (std::fabs(nume_b) < epsilon_parallel) )
81  {
82  return false; //COINCIDENT;
83  }
84  return false; //PARALLEL;
85  }
86 
87  double ua = nume_a / denom;
88  double ub = nume_b / denom;
89 
90  if( (ua >= 0.0) && (ua <= 1.0) && (ub >= 0.0) && (ub <= 1.0) )
91  {
92  /* // Get the intersection point. */
93  /* intersection[0] = second_segment_vertex[0][0] + */
94  /* ua*(second_segment_vertex[1][0] - second_segment_vertex[0][0]); */
95  /* intersection[1] = second_segment_vertex[0][1] + */
96  /* ua*(second_segment_vertex[1][1] - second_segment_vertex[0][1]); */
97 
98  return true; //INTERESECTING;
99  }
100 
101  return false; //NOT_INTERESECTING;
102  }
103 
104 }
105 
106 
110 
111 //=======================================================================
114 //=======================================================================
115 namespace SolarRadiationHelper
116  {
117 
118 
119 
120  //=======================================================================
122  //=======================================================================
123  void Zero_atmospheric_radiation_fct(const double& time,
124  double& solar_flux_magnitude,
125  Vector<double> &solar_flux_unit_vector,
126  double& total_diffuse_radiation)
127  {
128  solar_flux_magnitude=0.0;
129  solar_flux_unit_vector[0]= 0.0;
130  solar_flux_unit_vector[1]=-1.0;
131  total_diffuse_radiation=0.0;
132  }
133 
134 }
135 
139 
140 
141 
142 //======================================================================
145 //======================================================================
146 class SolarRadiationBase : public virtual FiniteElement,
147  public virtual FaceElement,
148  public virtual TemplateFreeUnsteadyHeatBaseFaceElement
149 
150 {
151 
152 public:
153 
155  SolarRadiationBase()
156  {
157  // Don't use tanh profile to smooth solar shadows
158  Smoothed_sun_shadow=false;
159 
160  // Factor for tanh profile to smooth solar shadows
161  Alpha_tanh_smooth_sun_shadow=100.0;
162 
163  // Assign default zero atmospheric radiation fct...
164  Atmospheric_radiation_fct_pt=
165  &SolarRadiationHelper::Zero_atmospheric_radiation_fct;
166 
167  // Make space for limiting angles for diffuse radiation
168  unsigned n_intpt = integral_pt()->nweight();
169  Diffuse_limit_angles.resize(n_intpt);
170  for (unsigned i=0;i<n_intpt;i++)
171  {
172  Diffuse_limit_angles[i].first=0.0;
173  Diffuse_limit_angles[i].second=MathematicalConstants::Pi;
174  }
175  }
176 
177 
180  void enable_smoothed_sun_shadow(const double& alpha_tanh_smooth_sun_shadow=
181  100.0)
182  {
183  Smoothed_sun_shadow=true;
184  Alpha_tanh_smooth_sun_shadow=alpha_tanh_smooth_sun_shadow;
185  }
186 
188  void disable_smoothed_sun_shadow()
189  {
190  Smoothed_sun_shadow=false;
191  }
192 
194  bool smoothed_sun_shadow_is_enabled()
195  {
196  return Smoothed_sun_shadow;
197  }
198 
199 
202  double alpha_tanh_smooth_sun_shadow()
203  {
204  return Alpha_tanh_smooth_sun_shadow;
205  }
206 
207 
209  void (* &atmospheric_radiation_fct_pt())(const double& time,
210  double &solar_flux_magnitude,
211  Vector<double>&
212  solar_flux_unit_vector,
213  double& total_diffuse_radiation)
214  {return Atmospheric_radiation_fct_pt;}
215 
216 
221  void update_limiting_angles(const Vector<Node*>&
222  shielding_node_pt);
223 
224 
228  virtual double atmospheric_radiation(const unsigned& intpt,
229  const double& time,
230  const Vector<double>& x,
231  const Vector<double>& n);
232 
233 
234 
236  void output_diffuse_radiation_cone(std::ostream &outfile,
237  const double& radius);
238 
239 
242 void output_diffuse_radiation_cone_min_angle(std::ostream &outfile,
243  const double& radius);
244 
247 void output_diffuse_radiation_cone_max_angle(std::ostream &outfile,
248  const double& radius);
249 
250 
252 void output_limiting_angles(std::ostream &outfile);
253 
254 
256 void output_atmospheric_radiation(std::ostream &outfile);
257 
258 protected:
259 
261  bool Smoothed_sun_shadow;
262 
264  double Alpha_tanh_smooth_sun_shadow;
265 
270  Vector<std::pair<double,double> > Diffuse_limit_angles;
271 
276  void (*Atmospheric_radiation_fct_pt)(const double& time,
277  double& solar_flux_magnitude,
278  Vector<double> &solar_flux_unit_vector,
279  double& total_diffuse_radiation);
280 
281  private:
282 
290  void check_quadrant_jump(const double& x_prev,
291  const double& y_prev,
292  const double& x_next,
293  const double& y_next,
294  bool& crossing_quadrants,
295  int& n_winding_increment,
296  std::string& error_string,
297  bool& no_problem);
298 
299 };
300 
301 
302 
303 
304 //=====================================================================
308 //=====================================================================
309 double SolarRadiationBase::atmospheric_radiation(const unsigned& intpt,
310  const double& time,
311  const Vector<double>& x,
312  const Vector<double>& n)
313 {
314  double radiation=0.0;
315  unsigned n_dim = this->nodal_dimension();
316  double solar_flux_magnitude=0.0;
317  Vector<double> solar_flux_unit_vector(n_dim);
318  double total_diffuse_radiation=0.0;
319  Atmospheric_radiation_fct_pt(time,solar_flux_magnitude,
320  solar_flux_unit_vector,
321  total_diffuse_radiation);
322 
323  // Diffuse radiation from opening angle
324  //-------------------------------------
325  double phi_exposed=Diffuse_limit_angles[intpt].second-
326  Diffuse_limit_angles[intpt].first;
327  if (phi_exposed>0.0)
328  {
329  radiation+=total_diffuse_radiation*phi_exposed/
330  MathematicalConstants::Pi;
331 
332 #ifdef PARANOID
333  double tol=1.0e-5;
334  if (phi_exposed-MathematicalConstants::Pi>tol)
335  {
336  std::stringstream error_message;
337  error_message << "Exposure angle " << phi_exposed
338  << " greater than 180^o at ( "
339  << x[0] << " , " << x[1] << " )\n"
340  << "Actual difference is: "
341  << phi_exposed-MathematicalConstants::Pi
342  << " which exceeds tolerance " << tol << std::endl;
343  //throw
344  OomphLibError(error_message.str(),
345  OOMPH_CURRENT_FUNCTION,
346  OOMPH_EXCEPTION_LOCATION);
347  }
348 #endif
349  }
350 
351  // Direct radiation (with directional cosine and ignore
352  //-----------------------------------------------------
353  // shielded bits)
354  //---------------
355 
356 
357  // Cos of angle between outer unit normal and direct solar flux
358  double cos_angle=-(n[0]*solar_flux_unit_vector[0]+
359  n[1]*solar_flux_unit_vector[1]);
360 
361  // No further contributions if we're pointing away from the sun
362  if (cos_angle>0.0)
363  {
364  // Angle that the incoming solar radiation makes
365  // with the horizontal
366  double theta=atan2(solar_flux_unit_vector[0],-solar_flux_unit_vector[1]);
367  double phi=0.5*MathematicalConstants::Pi+theta;
368 
369 
370  // Bounding angles of visibility
371  double phi_min=Diffuse_limit_angles[intpt].first;
372  double phi_max=Diffuse_limit_angles[intpt].second;
373 
374  // Is the sun visible?
375  if (Smoothed_sun_shadow)
376  {
377  // Smooth shadow
378  double tanh_factor=
379  0.5*(-tanh(Alpha_tanh_smooth_sun_shadow*(phi-phi_max))
380  -tanh(Alpha_tanh_smooth_sun_shadow*(phi_min-phi)));
381  radiation+=cos_angle*solar_flux_magnitude*tanh_factor;
382  }
383  else
384  {
385  // Hard shadow
386  if (phi>phi_min)
387  {
388  if (phi<phi_max)
389  {
390  radiation+=cos_angle*solar_flux_magnitude;
391  }
392  }
393  }
394  }
395 
396  return radiation;
397 }
398 
399 
400 //====================================================================
408 //====================================================================
409 void SolarRadiationBase::check_quadrant_jump(const double& x_prev,
410  const double& y_prev,
411  const double& x_next,
412  const double& y_next,
413  bool& crossing_quadrants,
414  int& n_winding_increment,
415  std::string& error_string,
416  bool& no_problem)
417 {
418  // Sanity check: We have to keep track of winding numbers
419  // but we're stuffed if we our discrete increments in coordinates
420  // jump diagonally across quadrants...
421  std::stringstream error_message;
422  no_problem=true;
423  crossing_quadrants=false;
424 
425  // y coordinate of intersection with x axis
426  double y_intersect=0.0;
427  double denom=x_next-x_prev;
428 
429  // If denominator is zero then we're intersection exactly
430  // at the sample point and we're absolutely screwed. This
431  // really shouldn't happen.
432  if (denom!=0.0)
433  {
434  y_intersect=y_prev-(y_next-y_prev)*x_prev/denom;
435  }
436  if ((x_prev>0.0)&&(y_prev>0.0)&&(x_next<0.0)&&(y_next<0.0))
437  {
438  crossing_quadrants=true;
439  error_message << "Jumped from upper right quadrant to lower left one\n";
440  oomph_info << "Jumped from upper right quadrant to lower left one\n";
441  if (y_intersect==0.0)
442  {
443  error_message << "..and cannot be repaired!\n";
444  no_problem=false;
445  }
446  // Path crosses the positive x axis: No winding
447  else if (y_intersect<0.0)
448  {
449  n_winding_increment=0;
450  }
451  // Path crosses the negative x axis from above: increase winding number
452  else if (y_intersect>0.0)
453  {
454  n_winding_increment=1;
455  }
456  }
457  else if ((x_prev<0.0)&&(y_prev<0.0)&&(x_next>0.0)&&(y_next>0.0))
458  {
459  crossing_quadrants=true;
460  error_message << "Jumped from lower left quadrant to upper right one\n";
461  oomph_info << "Jumped from lower left quadrant to upper right one\n";
462  if (y_intersect==0.0)
463  {
464  error_message << "..and cannot be repaired!\n";
465  no_problem=false;
466  }
467  // Path crosses the positive x axis: No winding
468  else if (y_intersect<0.0)
469  {
470  n_winding_increment=0;
471  }
472  // Path crosses the negative x axis from below: reduce winding number
473  else if (y_intersect>0.0)
474  {
475  n_winding_increment=-1;
476  }
477  }
478  else if ((x_prev<0.0)&&(y_prev>0.0)&&(x_next>0.0)&&(y_next<0.0))
479  {
480  crossing_quadrants=true;
481  error_message << "Jumped from upper left quadrant to lower right one\n";
482  oomph_info << "Jumped from upper left quadrant to lower right one\n";
483  if (y_intersect==0.0)
484  {
485  error_message << "..and cannot be repaired!\n";
486  no_problem=false;
487  }
488  // Path crosses the positive x axis: No winding
489  else if (y_intersect>0.0)
490  {
491  n_winding_increment=0;
492  }
493  // Path crosses the negative x axis from above: increase winding number
494  else if (y_intersect<0.0)
495  {
496  n_winding_increment=1;
497  }
498  }
499  else if ((x_prev>0.0)&&(y_prev<0.0)&&(x_next<0.0)&&(y_next>0.0))
500  {
501  crossing_quadrants=true;
502  error_message << "Jumped from lower right quadrant to upper left one\n";
503  oomph_info << "Jumped from lower right quadrant to upper left one\n";
504  if (y_intersect==0.0)
505  {
506  error_message << "..and cannot be repaired!\n";
507  no_problem=false;
508  }
509  // Path crosses the positive x axis: No winding
510  else if (y_intersect>0.0)
511  {
512  n_winding_increment=0;
513  }
514  // Path crosses the negative x axis from below: decrease winding number
515  else if (y_intersect<0.0)
516  {
517  n_winding_increment=-1;
518  }
519  }
520 }
521 
522 //=====================================================================
527 //=====================================================================
528  void SolarRadiationBase::update_limiting_angles(const Vector<Node*>&
529  shielding_node_pt)
530 {
531 
532  // Search through all shielding nodes and find the
533  // bounding ones for this element
534  Node* first_vertex_node_pt=node_pt(0);
535  unsigned nnod_el=nnode();
536  Node* second_vertex_node_pt=node_pt(nnod_el-1);
537 
538  // Find left and rightmost node of element in shielding_node_pt vector:
539  unsigned j_left=0;
540  bool found=false;
541  unsigned nnod=shielding_node_pt.size();
542  for (unsigned j=0;j<nnod;j++)
543  {
544  // We come from the left so we're definitely meeting the leftmost
545  // node
546  if ((shielding_node_pt[j]==first_vertex_node_pt)||
547  (shielding_node_pt[j]==second_vertex_node_pt))
548  {
549  j_left=j;
550  found=true;
551  break;
552  }
553  }
554  unsigned j_right=j_left+1;
555 
556  // Did we succeed?
557  if (!found)
558  {
559  throw OomphLibError("Failed to find leftmost node in shielding_node_pt",
560  OOMPH_CURRENT_FUNCTION,
561  OOMPH_EXCEPTION_LOCATION);
562  }
563 
564  //Set the value of n_intpt
565  unsigned n_intpt = integral_pt()->nweight();
566 
567  // Spatial dimension of the nodes
568  unsigned n_dim = this->nodal_dimension();
569  Vector<double> x(n_dim);
570  Vector<double> s(n_dim-1);
571 
572  //Loop over the integration points
573  for(unsigned ipt=0;ipt<n_intpt;ipt++)
574  {
575 
576  // Local coordinate of integration point
577  for (unsigned i=0;i<n_dim-1;i++)
578  {
579  s[i]=integral_pt()->knot(ipt,i);
580  }
581 
582  // No recycling of shape fcts -- may as well call interpolated_x
583  // directly
584  this->interpolated_x(s,x);
585 
586  oomph_info << "Integration point: " << x[0] << " " << x[1] << std::endl;
587 
588  // Loop over all potential shielding nodes to the left to find
589  // minimum angle
590  int n_winding=0;
591 
592  // Initial point
593  Node* nod_pt=shielding_node_pt[0];
594 
595  // Coordinate relative to sampling point
596  double x_prev=nod_pt->x(0)-x[0];
597  double y_prev=nod_pt->x(1)-x[1];
598 
599 #ifdef PARANOID
600  // Check that initial point is to the left of current point
601  // for calibration purposes...
602  if (x_prev>0.0)
603  {
604  throw OomphLibError(
605  "Leftmost point in shielding nodes is not to the left of current int pt",
606  OOMPH_CURRENT_FUNCTION,
607  OOMPH_EXCEPTION_LOCATION);
608  }
609 #endif
610 
611  // Angle against x-axis relative to sampling point
612  double phi_prev=atan2(y_prev,x_prev);
613 
614  oomph_info << "Vertex point " << 0 << " : "
615  << nod_pt->x(0) << " " << nod_pt->x(1)
616  << " phi_raw: " << phi_prev;
617  oomph_info << " wind1 " << n_winding << " ";
618 
619  // Want a positive angle!
620  if (phi_prev<0.0)
621  {
622  n_winding+=1;
623  phi_prev+=2.0*MathematicalConstants::Pi;
624  }
625 
626  oomph_info << " wind2 " << n_winding << " ";
627  oomph_info << " phi_adj: " << phi_prev << std::endl;
628 
629 
630  // Current best guess for minimum angle to the left
631  double phi_left=phi_prev;
632 
633 
634 
635  // Loop over all other nodes
636  for (unsigned j=1;j<=j_left;j++)
637  {
638  Node* nod_pt=shielding_node_pt[j];
639  double x_next=nod_pt->x(0)-x[0];
640  double y_next=nod_pt->x(1)-x[1];
641  double phi_next=atan2(y_next,x_next);
642 
643  oomph_info << "Vertex point " << j << " : "
644  << nod_pt->x(0) << " " << nod_pt->x(1)
645  << " phi_raw: " << phi_next;
646 
647  // Sanity check: We have to keep track of winding numbers
648  // but we're stuffed if we our discrete increments in coordinates
649  // jump diagonally across quadrants...
650  std::string error_string;
651  bool no_problem=true;
652  bool crossing_quadrants=false;
653  int n_winding_increment=0;
654 
655  // Check if we've crossed quadrants
656  check_quadrant_jump(x_prev,
657  y_prev,
658  x_next,
659  y_next,
660  crossing_quadrants,
661  n_winding_increment,
662  error_string,
663  no_problem);
664 
665  // Success?
666  if (no_problem)
667  {
668  n_winding+=n_winding_increment;
669  }
670  // Complain bitterly
671  else
672  {
673  std::stringstream error_message;
674  error_message << error_string;
675  error_message << "\n x/y_prev, x/y_next:"
676  << x_prev << " "
677  << y_prev << " "
678  << x_next << " "
679  << y_next << " "<< std::endl;
680  error_message << "for point at "
681  << x[0] << " " << x[1] << std::endl;
682  error_message << "chain of shielding nodes on left:\n";
683  for (unsigned jj=1;jj<=j_left;jj++)
684  {
685  Node* nnod_pt=shielding_node_pt[jj];
686  error_message << nnod_pt->x(0) << " "
687  << nnod_pt->x(1) << "\n";
688  }
689  throw OomphLibError(error_message.str(),
690  OOMPH_CURRENT_FUNCTION,
691  OOMPH_EXCEPTION_LOCATION);
692  }
693 
694  oomph_info << " wind1 " << n_winding << " ";
695 
696  // If we're crossing the x axis to the left of the origin
697  // without jumping across quadrants we
698  // have to adjust the winding number
699  if (!crossing_quadrants)
700  {
701  if ((x_prev<0.0)&&(x_next<0.0))
702  {
703  if ((y_prev>=0.0)&&(y_next<0.0))
704  {
705  n_winding+=1;
706  }
707  else if ((y_prev<0.0)&&(y_next>=0.0))
708  {
709  n_winding-=1;
710  }
711  }
712  }
713 
714 
715  // Correct angle by winding
716  phi_next+=double(n_winding)*2.0*MathematicalConstants::Pi;
717 
718  oomph_info << " wind2 " << n_winding << " ";
719  oomph_info << " phi_adj: " << phi_next << std::endl;
720 
721  // Is it smaller?
722  if (phi_next<phi_left)
723  {
724  phi_left=phi_next;
725  }
726 
727  // Bump up
728  x_prev=x_next;
729  y_prev=y_next;
730  phi_prev=phi_next;
731  }
732 
733 
734  // Loop over all potential shielding nodes to the right to find
735  // maximum angle
736  n_winding=0;
737 
738  // Initial point
739  nod_pt=shielding_node_pt[nnod-1];
740 
741  // Coordinate relative to sampling point
742  x_prev=nod_pt->x(0)-x[0];
743  y_prev=nod_pt->x(1)-x[1];
744 
745 #ifdef PARANOID
746  // Check that initial point is to the right of current point
747  // for calibration purposes...
748  if (x_prev<0.0)
749  {
750  throw OomphLibError(
751  "Rightmost point in shielding nodes is not to the right of current int pt",
752  OOMPH_CURRENT_FUNCTION,
753  OOMPH_EXCEPTION_LOCATION);
754  }
755 #endif
756 
757  // Angle against x-axis relative to sampling point
758  phi_prev=atan2(y_prev,x_prev);
759 
760  // Current best guess for maximum angle to the right
761  double phi_right=phi_prev;
762 
763  // Loop over all other nodes
764  //for (unsigned j=j_right+1;j<nnod;j++)
765  for (unsigned j=nnod-2;j>=j_right;j--)
766  {
767  Node* nod_pt=shielding_node_pt[j];
768  double x_next=nod_pt->x(0)-x[0];
769  double y_next=nod_pt->x(1)-x[1];
770  double phi_next=atan2(y_next,x_next);
771 
772  // Sanity check: We have to keep track of winding numbers
773  // but we're stuffed if we our discrete increments in coordinates
774  // jump diagonally across quadrants...
775  std::string error_string;
776  bool no_problem=true;
777  bool crossing_quadrants=false;
778  int n_winding_increment=0;
779 
780  // Check if we've crossed quadrants
781  check_quadrant_jump(x_prev,
782  y_prev,
783  x_next,
784  y_next,
785  crossing_quadrants,
786  n_winding_increment,
787  error_string,
788  no_problem);
789 
790  // Success?
791  if (no_problem)
792  {
793  n_winding+=n_winding_increment;
794  }
795  // Complain bitterly
796  else
797  {
798  std::stringstream error_message;
799  error_message << error_string;
800  error_message << "\n x/y_prev, x/y_next:"
801  << x_prev << " "
802  << y_prev << " "
803  << x_next << " "
804  << y_next << " "<< std::endl;
805  error_message << "for point at "
806  << x[0] << " " << x[1] << std::endl;
807  error_message << "chain of shielding nodes on right:\n";
808  for (unsigned jj=j_right;jj<nnod;jj++)
809  {
810  Node* nnod_pt=shielding_node_pt[jj];
811  error_message << nnod_pt->x(0) << " "
812  << nnod_pt->x(1) << "\n";
813  }
814  throw OomphLibError(error_message.str(),
815  OOMPH_CURRENT_FUNCTION,
816  OOMPH_EXCEPTION_LOCATION);
817  }
818 
819 
820  // If we're crossing the x axis to the left of the origin
821  // without jumping across quadrants we
822  // have to adjust the winding number
823  if (!crossing_quadrants)
824  {
825  if ((x_prev<0.0)&&(x_next<0.0))
826  {
827  if ((y_prev>=0.0)&&(y_next<0.0))
828  {
829  n_winding+=1;
830  }
831  else if ((y_prev<0.0)&&(y_next>=0.0))
832  {
833  n_winding-=1;
834  }
835  }
836  }
837 
838  // Correct angle by winding
839  phi_next+=double(n_winding)*2.0*MathematicalConstants::Pi;
840 
841  // Is it bigger?
842  if (phi_next>phi_right)
843  {
844  phi_right=phi_next;
845  }
846 
847  // Bump up
848  x_prev=x_next;
849  y_prev=y_next;
850  phi_prev=phi_next;
851  }
852 
853 //---
854 
855  Diffuse_limit_angles[ipt].first=phi_right;
856  Diffuse_limit_angles[ipt].second=phi_left;
857  }
858 }
859 
860 
861 
862 //=====================================================================
864 //=====================================================================
865 void SolarRadiationBase::output_diffuse_radiation_cone(
866  std::ostream &outfile, const double& radius)
867 {
868  //Set the value of n_intpt
869  unsigned n_intpt = integral_pt()->nweight();
870 
871  // Spatial dimension of the nodes
872  unsigned n_dim = this->nodal_dimension();
873  Vector<double> x(n_dim);
874  Vector<double> s(n_dim-1);
875 
876  //Loop over the integration points
877  for(unsigned ipt=0;ipt<n_intpt;ipt++)
878  {
879  outfile << "ZONE\n";
880 
881  // Local coordinate of integration point
882  for (unsigned i=0;i<n_dim-1;i++)
883  {
884  s[i]=integral_pt()->knot(ipt,i);
885  }
886 
887  // No recycling of shape fcts -- may as well call interpolated_x
888  // directly
889  this->interpolated_x(s,x);
890 
891  double phi_min=Diffuse_limit_angles[ipt].first;
892  double phi_max=Diffuse_limit_angles[ipt].second;
893 
894  if (phi_max>phi_min)
895  {
896  outfile << x[0]+radius*cos(phi_min) << " "
897  << x[1]+radius*sin(phi_min) << "\n"
898  << x[0] << " "
899  << x[1] << "\n"
900  << x[0]+radius*cos(phi_max) << " "
901  << x[1]+radius*sin(phi_max) << "\n";
902  }
903  else
904  {
905  outfile << x[0] << " "
906  << x[1] << "\n"
907  << x[0] << " "
908  << x[1] << "\n"
909  << x[0] << " "
910  << x[1] << "\n";
911  }
912  }
913 }
914 
915 
916 //=====================================================================
919 //=====================================================================
920 void SolarRadiationBase::output_diffuse_radiation_cone_max_angle(
921  std::ostream &outfile, const double& radius)
922 {
923  //Set the value of n_intpt
924  unsigned n_intpt = integral_pt()->nweight();
925 
926  // Spatial dimension of the nodes
927  unsigned n_dim = this->nodal_dimension();
928  Vector<double> x(n_dim);
929  Vector<double> s(n_dim-1);
930 
931  //Loop over the integration points
932  for(unsigned ipt=0;ipt<n_intpt;ipt++)
933  {
934  outfile << "ZONE\n";
935 
936  // Local coordinate of integration point
937  for (unsigned i=0;i<n_dim-1;i++)
938  {
939  s[i]=integral_pt()->knot(ipt,i);
940  }
941 
942  // No recycling of shape fcts -- may as well call interpolated_x
943  // directly
944  this->interpolated_x(s,x);
945 
946  double phi_max=Diffuse_limit_angles[ipt].second;
947 
948  unsigned nplot=100;
949  double fract=0.1;
950  for (unsigned j=0;j<nplot;j++)
951  {
952  double phi=phi_max*double(j)/double(nplot-1);
953 
954  outfile
955  << x[0]+fract*(1.0+0.1*double(j)/double(nplot-1))*radius*cos(phi) << " "
956  << x[1]+fract*(1.0+0.1*double(j)/double(nplot-1))*radius*sin(phi) << "\n";
957  }
958  outfile << x[0] << " "
959  << x[1] << "\n"
960  << x[0]+radius*cos(phi_max) << " "
961  << x[1]+radius*sin(phi_max) << "\n";
962 
963  }
964 }
965 //=====================================================================
968 //=====================================================================
969 void SolarRadiationBase::output_diffuse_radiation_cone_min_angle(
970  std::ostream &outfile, const double& radius)
971 {
972  //Set the value of n_intpt
973  unsigned n_intpt = integral_pt()->nweight();
974 
975  // Spatial dimension of the nodes
976  unsigned n_dim = this->nodal_dimension();
977  Vector<double> x(n_dim);
978  Vector<double> s(n_dim-1);
979 
980  //Loop over the integration points
981  for(unsigned ipt=0;ipt<n_intpt;ipt++)
982  {
983  outfile << "ZONE\n";
984 
985  // Local coordinate of integration point
986  for (unsigned i=0;i<n_dim-1;i++)
987  {
988  s[i]=integral_pt()->knot(ipt,i);
989  }
990 
991  // No recycling of shape fcts -- may as well call interpolated_x
992  // directly
993  this->interpolated_x(s,x);
994 
995  double phi_min=Diffuse_limit_angles[ipt].first;
996 
997  unsigned nplot=100;
998  double fract=0.13;
999  for (unsigned j=0;j<nplot;j++)
1000  {
1001  double phi=phi_min*double(j)/double(nplot-1);
1002 
1003  outfile
1004  << x[0]+fract*(1.0+0.1*double(j)/double(nplot-1))*radius*cos(phi) << " "
1005  << x[1]+fract*(1.0+0.1*double(j)/double(nplot-1))*radius*sin(phi) << "\n";
1006  }
1007  outfile << x[0] << " "
1008  << x[1] << "\n"
1009  << x[0]+radius*cos(phi_min) << " "
1010  << x[1]+radius*sin(phi_min) << "\n";
1011 
1012  }
1013 }
1014 
1015 
1016 
1017 //=====================================================================
1019 //=====================================================================
1020 void SolarRadiationBase::output_limiting_angles(std::ostream &outfile)
1021 {
1022  //Set the value of n_intpt
1023  unsigned n_intpt = integral_pt()->nweight();
1024 
1025  // Spatial dimension of the nodes
1026  unsigned n_dim = this->nodal_dimension();
1027  Vector<double> x(n_dim);
1028  Vector<double> s(n_dim-1);
1029  Vector<double> normal(n_dim);
1030 
1031  //Loop over the integration points
1032  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1033  {
1034  outfile << "ZONE\n";
1035 
1036  // Local coordinate of integration point
1037  for (unsigned i=0;i<n_dim-1;i++)
1038  {
1039  s[i]=integral_pt()->knot(ipt,i);
1040  }
1041 
1042  // No recycling of shape fcts -- may as well call interpolated_x
1043  // directly
1044  this->interpolated_x(s,x);
1045 
1046  // Get outer unit normal
1047  this->outer_unit_normal(s,normal);
1048 
1049  double phi_min=Diffuse_limit_angles[ipt].first;
1050  double phi_max=Diffuse_limit_angles[ipt].second;
1051 
1052  outfile << x[0] << " "
1053  << x[1] << " "
1054  << normal[0] << " "
1055  << normal[1] << " "
1056  << phi_min << " "
1057  << phi_max << "\n";
1058 
1059  }
1060 }
1061 
1062 
1063 
1064 
1065 //=====================================================================
1067 //=====================================================================
1068 void SolarRadiationBase::output_atmospheric_radiation(std::ostream &outfile)
1069 {
1070 
1071  // Get time from first node
1072  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1073 
1074  //Set the value of n_intpt
1075  unsigned n_intpt = integral_pt()->nweight();
1076 
1077  // Spatial dimension of the nodes
1078  unsigned n_dim = this->nodal_dimension();
1079  Vector<double> x(n_dim);
1080  Vector<double> s(n_dim-1);
1081  Vector<double> normal(n_dim);
1082 
1083  //Loop over the integration points
1084  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1085  {
1086  outfile << "ZONE\n";
1087 
1088  // Local coordinate of integration point
1089  for (unsigned i=0;i<n_dim-1;i++)
1090  {
1091  s[i]=integral_pt()->knot(ipt,i);
1092  }
1093 
1094  // No recycling of shape fcts -- may as well call interpolated_x
1095  // directly
1096  this->interpolated_x(s,x);
1097 
1098  // Get outer unit normal
1099  this->outer_unit_normal(s,normal);
1100 
1101  // Get radiation
1102  double radiation=atmospheric_radiation(ipt,
1103  time,
1104  x,
1105  normal);
1106  // output
1107  outfile << x[0] << " "
1108  << x[1] << " "
1109  << radiation << " "
1110  << normal[0] << " "
1111  << normal[1] << "\n";
1112 
1113  }
1114 }
1115 
1116 
1117 
1121 
1122 
1123 
1124 //======================================================================
1127 //======================================================================
1128 class StefanBoltzmannRadiationBase :
1129 public virtual FiniteElement,
1130  public virtual FaceElement,
1131  public virtual TemplateFreeUnsteadyHeatBaseFaceElement
1132 {
1133 
1134 public:
1135 
1137  StefanBoltzmannRadiationBase()
1138  {
1139  // Pointer to non-dim Stefan Boltzmann constant
1140  Sigma_pt=0;
1141 
1142  // Pointer to non-dim zero centrigrade offset in Stefan Boltzmann law
1143  Theta_0_pt=0;
1144 
1145  // Make space for Stefan Boltzmann illumination info
1146  unsigned n_intpt = integral_pt()->nweight();
1147  Stefan_boltzmann_illumination_info.resize(n_intpt);
1148  }
1149 
1151  double*& sigma_pt(){return Sigma_pt;}
1152 
1154  double*& theta_0_pt(){return Theta_0_pt;}
1155 
1157  double sigma()
1158  {
1159  if (Sigma_pt==0)
1160  {
1161  return 0.0;
1162  }
1163  else
1164  {
1165  return *Sigma_pt;
1166  }
1167  }
1168 
1172  double theta_0()
1173  {
1174  if (Sigma_pt==0)
1175  {
1176  return 0.0;
1177  }
1178  else
1179  {
1180  if (Theta_0_pt==0)
1181  {
1182  throw OomphLibError("Non-dim zero-centrigrade offset not set",
1183  OOMPH_CURRENT_FUNCTION,
1184  OOMPH_EXCEPTION_LOCATION);
1185  }
1186  else
1187  {
1188  return *Theta_0_pt;
1189  }
1190  }
1191  }
1192 
1193 
1199  double contribution_to_stefan_boltzmann_radiation(
1200  const Vector<double>& r_illuminated,
1201  const Vector<double>& n_illuminated,
1202  const Vector<unsigned>& visible_intpts_in_current_element)
1203  {
1204  // Dummy file
1205  std::ofstream outfile;
1206  return contribution_to_stefan_boltzmann_radiation(
1207  r_illuminated,
1208  n_illuminated,
1209  visible_intpts_in_current_element,
1210  outfile);
1211  }
1212 
1218  double contribution_to_stefan_boltzmann_radiation(
1219  const Vector<double>& r_illuminated,
1220  const Vector<double>& n_illuminated,
1221  const Vector<unsigned>& visible_intpts_in_current_element,
1222  std::ofstream& outfile);
1223 
1225  void wipe_stefan_boltzmann_illumination_info()
1226  {
1227  unsigned nint=Stefan_boltzmann_illumination_info.size();
1228  for (unsigned ipt=0;ipt<nint;ipt++)
1229  {
1230  Stefan_boltzmann_illumination_info[ipt].clear();
1231  }
1232  }
1233 
1247  void add_stefan_boltzmann_illumination_info(
1248  const unsigned& ipt,
1249  StefanBoltzmannRadiationBase* illuminating_el_pt,
1250  Vector<unsigned>& illuminating_integration_point_index,
1251  const bool& add_solid_position_data=true)
1252  {
1253 #ifdef PARANOID
1254  unsigned n=Stefan_boltzmann_illumination_info[ipt].size();
1255  for (unsigned e=0;e<n;e++)
1256  {
1257  if (Stefan_boltzmann_illumination_info[ipt][e].first==illuminating_el_pt)
1258  {
1259  throw OomphLibError(
1260  "Element has already been added!",
1261  OOMPH_CURRENT_FUNCTION,
1262  OOMPH_EXCEPTION_LOCATION);
1263  }
1264  }
1265 #endif
1266 
1267  Stefan_boltzmann_illumination_info[ipt].push_back(
1268  std::make_pair(illuminating_el_pt,illuminating_integration_point_index));
1269 
1270  // Add nodal Data of illuminating elements as external data
1271  unsigned nnod=illuminating_el_pt->nnode();
1272  for (unsigned j=0;j<nnod;j++)
1273  {
1274  Node* ext_node_pt=illuminating_el_pt->node_pt(j);
1275  bool own=false;
1276  for (unsigned jj=0;jj<nnod;jj++)
1277  {
1278  if (node_pt(jj)==ext_node_pt)
1279  {
1280  own=true;
1281  break;
1282  }
1283  }
1284  if (!own)
1285  {
1286  add_external_data(ext_node_pt);
1287  if (add_solid_position_data)
1288  {
1289  SolidNode* solid_node_pt=dynamic_cast<SolidNode*>(ext_node_pt);
1290  if (solid_node_pt!=0)
1291  {
1292  add_external_data(solid_node_pt->variable_position_pt());
1293  }
1294  }
1295  }
1296  }
1297 
1298  }
1299 
1301  void output_stefan_boltzmann_radiation(std::ostream &outfile);
1302 
1305  void output_stefan_boltzmann_radiation_rays(
1306  std::ostream &outfile,
1307  const unsigned& integration_point=UINT_MAX);
1308 
1309 
1314  double incoming_stefan_boltzmann_radiation(const unsigned& ipt,
1315  const Vector<double>& r_illuminated,
1316  const Vector<double>& n_illuminated)
1317  {
1318  // Initialise flux
1319  double flux=0.0;
1320 
1322  unsigned n_contrib=Stefan_boltzmann_illumination_info[ipt].size();
1323  for (unsigned e=0;e<n_contrib;e++)
1324  {
1325  // Pointer to contributing (illuminating) element
1326  StefanBoltzmannRadiationBase* el_pt=
1327  Stefan_boltzmann_illumination_info[ipt][e].first;
1328 
1329  // Indices of illuminating integration points that contribute
1330  Vector<unsigned> visible_integration_points=
1331  Stefan_boltzmann_illumination_info[ipt][e].second;
1332 
1333  // Add contribution
1334  flux+=el_pt->contribution_to_stefan_boltzmann_radiation
1335  (r_illuminated,n_illuminated,visible_integration_points);
1336  }
1337 
1338  return flux;
1339  }
1340 
1341  protected:
1342 
1344  double* Sigma_pt;
1345 
1347  double* Theta_0_pt;
1348 
1360  Vector<Vector<std::pair<StefanBoltzmannRadiationBase*,Vector<unsigned> > > >
1361  Stefan_boltzmann_illumination_info;
1362 
1363 
1364 };
1365 
1366 
1367 
1368 //=====================================================================
1373 //=====================================================================
1374 double StefanBoltzmannRadiationBase::contribution_to_stefan_boltzmann_radiation(
1375  const Vector<double>& r_illuminated,
1376  const Vector<double>& n_illuminated,
1377  const Vector<unsigned>& visible_intpts_in_current_element,
1378  std::ofstream &outfile)
1379 {
1380 
1381  if (outfile.is_open())
1382  {
1383  outfile << "ZONE\n";
1384  outfile << r_illuminated[0] << " "
1385  << r_illuminated[1] << "0 0\n";
1386  }
1387 
1388  // Initialise
1389  double contribution=0.0;
1390 
1391  //Find out how many nodes there are
1392  unsigned n_node = nnode();
1393 
1394  //Set up memory for the shape functions
1395  Shape psi(n_node);
1396  DShape dpsids(n_node,1);
1397 
1398  // Loop over contributing integration points
1399  unsigned nint=visible_intpts_in_current_element.size();
1400  for (unsigned ii=0;ii<nint;ii++)
1401  {
1402  // Get integration point
1403  unsigned ipt=visible_intpts_in_current_element[ii];
1404 
1405  //Only need to call the local derivatives
1406  dshape_local_at_knot(ipt,psi,dpsids);
1407 
1408  //Get the integral weight
1409  double w = integral_pt()->weight(ipt);
1410 
1411  //Calculate the coords and tangent vector
1412  Vector<double> r_illuminating(2,0.0);
1413  Vector<double> interpolated_dxds(2,0.0);
1414  double interpolated_u=0;
1415 
1416  //Loop over nodes
1417  for(unsigned l=0;l<n_node;l++)
1418  {
1419  // Add to temperature
1420  interpolated_u+=
1421  node_pt(l)->value(U_index_ust_heat)*psi[l];
1422 
1423  //Loop over directions
1424  for(unsigned i=0;i<2;i++)
1425  {
1426  //ALH: Commented out because this is the memory error
1427  //s[i]=this->integral_pt()->knot(ipt,i);
1428  r_illuminating[i] += nodal_position(l,i)*psi(l);
1429  interpolated_dxds[i] += nodal_position(l,i)*dpsids(l,0);
1430  }
1431  }
1432 
1433  // Vector from illuminated to illuminated point
1434  Vector<double> ray(2);
1435  ray[0]=r_illuminating[0]-r_illuminated[0];
1436  ray[1]=r_illuminating[1]-r_illuminated[1];
1437 
1438 
1439  // Get inverse length
1440  double inv_length=1.0/sqrt(ray[0]*ray[0]+ray[1]*ray[1]);
1441 
1442  // e_phi (phi measured in mathematically negative sense
1443  // from vertically above illuminated point -- 'cos that's
1444  // what it is in my sketch..
1445  Vector<double> e_phi(2);
1446  e_phi[0]= inv_length*ray[1];
1447  e_phi[1]=-inv_length*ray[0];
1448 
1449  // Angles
1450  double cos_phi=inv_length*(n_illuminated[0]*ray[0]+
1451  n_illuminated[1]*ray[1]);
1452 
1453  // INTEGRAL IS:
1454  // \int sigma T^4 1/2 \cos \varphi d \varphi =
1455  // where
1456  // d\varphi=1/|{\bf R}| \partial {\bf R}/\partial s \cdot {\bf e}_\varphi ds
1457 
1458  double integrand=
1459  this->sigma()*pow((interpolated_u+this->theta_0()),4)*
1460  0.5*cos_phi*std::fabs(inv_length*(interpolated_dxds[0]*e_phi[0]+
1461  interpolated_dxds[1]*e_phi[1]));
1462 
1463 
1464  // Add it in...
1465  contribution+=integrand*w;
1466 
1467  if (outfile.is_open())
1468  {
1469  outfile << r_illuminating[0] << " "
1470  << r_illuminating[1] << " "
1471  << cos_phi << " "
1472  << integrand << " "
1473  << "\n";
1474  }
1475  }
1476 
1477  return contribution;
1478 }
1479 
1480 
1481 
1482 
1483 //=====================================================================
1486 //=====================================================================
1487 void StefanBoltzmannRadiationBase::output_stefan_boltzmann_radiation_rays(
1488  std::ostream &outfile, const unsigned& integration_point)
1489 {
1490 
1491  // Vector to illuminated/ing Gauss points
1492  Vector<double> r_illuminated(2);
1493  Vector<double> s_illuminated(1);
1494  Vector<double> unit_normal_illuminated(2);
1495  Vector<double> r_illuminating(2);
1496  Vector<double> s_illuminating(1);
1497  Vector<double> unit_normal_illuminating(2);
1498 
1499  unsigned ipt_lo=0;
1500  unsigned ipt_hi=this->integral_pt()->nweight()-1;
1501  if (integration_point!=UINT_MAX)
1502  {
1503  ipt_lo=integration_point;
1504  ipt_hi=integration_point;
1505  }
1506 
1507  // Loop over integration points
1508  for (unsigned ipt=ipt_lo;ipt<=ipt_hi;ipt++)
1509  {
1510  // Local coordinate of integration point
1511  s_illuminated[0]=this->integral_pt()->knot(ipt,0);
1512 
1513  // Get coordinate of illuminated integration point
1514  this->interpolated_x(s_illuminated,r_illuminated);
1515 
1516  // Plot illuminted point
1517  outfile << "ZONE\n";
1518  outfile << r_illuminated[0] << " " << r_illuminated[1] << "\n";
1519 
1520  // Number of illuminating elements
1521  unsigned n_illuminating=Stefan_boltzmann_illumination_info[ipt].size();
1522  for (unsigned i2=0;i2<n_illuminating;i2++)
1523  {
1524  // Get pointer to illuminating element
1525  FiniteElement* el_pt=
1526  (//dynamic_cast<FiniteElement*>(
1527  Stefan_boltzmann_illumination_info[ipt][i2].first);
1528 
1529  // Illuminating Gauss points in illuminating element
1530  unsigned n_pts=(Stefan_boltzmann_illumination_info[ipt][i2].second).size();
1531  for (unsigned ipt2=0;ipt2<n_pts;ipt2++)
1532  {
1533  // Illuminating integration point
1534  unsigned ipt_illum=
1535  (Stefan_boltzmann_illumination_info[ipt][i2].second)[ipt2];
1536 
1537  // Get local coordinate of that integration point
1538  s_illuminating[0]=el_pt->integral_pt()->knot(ipt_illum,0);
1539 
1540  // Get coordinate of illuminating integration point
1541  el_pt->interpolated_x(s_illuminating,r_illuminating);
1542 
1543  // Plot ray to illuminting point and back
1544  outfile << r_illuminating[0] << " " << r_illuminating[1] << "\n";
1545  outfile << r_illuminated[0] << " " << r_illuminated[1] << "\n";
1546  }
1547  }
1548  }
1549 }
1550 
1551 
1552 
1556 
1557 
1558 
1559 //======================================================================
1564 //======================================================================
1565 template <class ELEMENT>
1566 class StefanBoltzmannUnsteadyHeatFluxElement :
1567 public virtual FaceGeometry<ELEMENT>,
1568  public virtual FaceElement,
1569  public virtual StefanBoltzmannRadiationBase,
1570  public virtual SolarRadiationBase,
1571  public virtual UnsteadyHeatBaseFaceElement<ELEMENT>
1572 {
1573 
1574 public:
1575 
1578  StefanBoltzmannUnsteadyHeatFluxElement(FiniteElement* const &bulk_el_pt,
1579  const int &face_index);
1580 
1582  StefanBoltzmannUnsteadyHeatFluxElement(
1583  const StefanBoltzmannUnsteadyHeatFluxElement& dummy)
1584  {
1585  BrokenCopy::broken_copy("StefanBoltzmannUnsteadyHeatFluxElement");
1586  }
1587 
1590  void set_integration_scheme(Integral* const &integral_pt)
1591  {
1592  FiniteElement::set_integration_scheme(integral_pt);
1593  unsigned n_intpt = integral_pt->nweight();
1594  Stefan_boltzmann_illumination_info.resize(n_intpt);
1595  }
1596 
1602  double zeta_nodal(const unsigned &n, const unsigned &k,
1603  const unsigned &i) const
1604  {return FaceElement::zeta_nodal(n,k,i);}
1605 
1606 
1613  void output_stefan_boltzmann_radiation(std::ostream &outfile)
1614  {
1615  unsigned n_dim = this->nodal_dimension();
1616  Vector<double> x(n_dim);
1617  Vector<double> s(n_dim-1);
1618  Vector<double> unit_normal(n_dim);
1619 
1620  // Get continuous time from timestepper of first node
1621  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1622 
1623  // Number of integration points
1624  unsigned n_intpt = integral_pt()->nweight();
1625 
1626  outfile << "ZONE\n";
1627 
1628  //Loop over the integration points
1629  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1630  {
1631  // Local coordinate of integration point
1632  for (unsigned i=0;i<n_dim-1;i++)
1633  {
1634  s[i]=integral_pt()->knot(ipt,i);
1635  }
1636 
1637  // Get Eulerian coordinates
1638  this->interpolated_x(s,x);
1639 
1640  // Get temperature
1641  double interpolated_u=0;
1642  this->interpolated_u(s,interpolated_u);
1643 
1644  // Outer unit normal
1645  outer_unit_normal(s,unit_normal);
1646 
1647  // Outgoing Stefan Boltzmann radiation
1648  double outgoing_sb_radiation=
1649  this->sigma()*pow((interpolated_u+this->theta_0()),4);
1650 
1651  // Incoming Stefan Boltzmann radiation
1652  double incoming_sb_radiation=incoming_stefan_boltzmann_radiation
1653  (ipt,x,unit_normal);
1654 
1655  // Get the net incoming flux (should be the difference between
1656  // the two previous ones!
1657  double flux=0.0;
1658  get_flux(ipt,time,x,unit_normal,interpolated_u,flux);
1659 
1660 
1661 // hierher paranoidify
1662 //#ifdef PARANOID
1663  if (fabs(flux-(incoming_sb_radiation-outgoing_sb_radiation)))
1664  {
1665  std::stringstream error_message;
1666  error_message
1667  << "Difference between incoming ( " << incoming_sb_radiation
1668  << " ) and outgoing ( " << outgoing_sb_radiation << " ) is "
1669  << (incoming_sb_radiation-outgoing_sb_radiation)
1670  << " but total flux is " << flux << " so difference is "
1671  << fabs(flux-(incoming_sb_radiation-outgoing_sb_radiation));
1672  throw OomphLibError(
1673  error_message.str(),
1674  OOMPH_CURRENT_FUNCTION,
1675  OOMPH_EXCEPTION_LOCATION);
1676  }
1677 //#endif
1678 
1679 
1680  //Output the x,y,..
1681  for(unsigned i=0;i<n_dim;i++)
1682  {
1683  outfile << x[i] << " ";
1684  }
1685 
1686  // Temperature
1687  outfile << interpolated_u << " ";
1688 
1689  // Total incoming flux (may include additional contributions)
1690  outfile << flux << " ";
1691 
1692  // Incoming sb flux
1693  outfile << incoming_sb_radiation << " ";
1694 
1695  // Outgoing sb flux
1696  outfile << outgoing_sb_radiation << " ";
1697 
1698  // Outer unit normal
1699  for(unsigned i=0;i<n_dim;i++)
1700  {
1701  outfile << unit_normal[i] << " ";
1702  }
1703  outfile << std::endl;
1704  }
1705  }
1706 
1707 
1708 
1709  protected:
1710 
1711 
1721  virtual void get_flux(const unsigned& ipt,
1722  const double& time,
1723  const Vector<double>& x,
1724  const Vector<double>& n,
1725  const double& u,
1726  double& flux)
1727  {
1728 
1729 #ifdef PARANOID
1730  // We shouldn't really have an externally imposed flux (cos it's ignored)
1731  if (Flux_fct_pt != 0)
1732  {
1733  //Issue a warning
1734  OomphLibWarning("Flux_fct_pt is ignored in this element\n",
1735  OOMPH_CURRENT_FUNCTION,
1736  OOMPH_EXCEPTION_LOCATION);
1737  }
1738 #endif
1739 
1740  // Outgoing Stefan Boltzmann radiation in terms of pre-computed
1741  // local temperature
1742  double outgoing_sb_radiation=this->sigma()*pow((u+this->theta_0()),4);
1743 
1744  // Incoming Stefan Boltzmann radiation
1745  double incoming_sb_radiation=incoming_stefan_boltzmann_radiation(ipt,x,n);
1746 
1747  // Atmospheric radiation
1748  double incoming_atmospheric=atmospheric_radiation(ipt,time,x,n);
1749 
1750  // Net flux
1751  flux=incoming_sb_radiation-outgoing_sb_radiation+incoming_atmospheric;
1752 
1753 
1754  }
1755 
1756 
1757 
1758 };
1759 
1763 
1764 //===========================================================================
1767 //===========================================================================
1768 template<class ELEMENT>
1769 StefanBoltzmannUnsteadyHeatFluxElement<ELEMENT>::
1770 StefanBoltzmannUnsteadyHeatFluxElement(FiniteElement* const &bulk_el_pt,
1771  const int &face_index) :
1772 UnsteadyHeatBaseFaceElement<ELEMENT>(bulk_el_pt,face_index)
1773 {
1774 
1775 #ifdef PARANOID
1776  {
1777  //Check that the element is not a refineable 3d element
1778  if(bulk_el_pt->dim()==3)
1779  {
1780  //Is it refineable
1781  if(dynamic_cast<RefineableElement*>(bulk_el_pt))
1782  {
1783  //Issue a warning
1784  std::string error_string=
1785  "This face element will not work correctly if nodes are hanging.\n";
1786  error_string+="Use the refineable version instead. ";
1787  OomphLibWarning(error_string,
1788  OOMPH_CURRENT_FUNCTION,
1789  OOMPH_EXCEPTION_LOCATION);
1790  }
1791  }
1792  }
1793 #endif
1794 
1795 #ifdef PARANOID
1796  // Check spatial dimension
1797  if (Dim!=2)
1798  {
1799  //Issue a warning
1800  throw OomphLibError(
1801  "This element will almost certainly not work in non-2D problems, though it should be easy enough to upgrade... Volunteers?\n",
1802  OOMPH_CURRENT_FUNCTION,
1803  OOMPH_EXCEPTION_LOCATION);
1804  }
1805 #endif
1806 
1807 }
1808 
1812 
1813 
1814 //=======================================================================
1816 //=======================================================================
1817 namespace StefanBoltzmannHelper
1818 {
1819 
1821  Vector<Vector<std::set<FiniteElement*> > > Element_in_bin;
1822 
1824  Vector<double> Max_coord(2,-DBL_MAX);
1825 
1827  Vector<double> Min_coord(2, DBL_MAX);
1828 
1830  Vector<double> Dx(2,0.0);
1831 
1833  unsigned Nx_bin=1000;
1834 
1836  unsigned Ny_bin=1000;
1837 
1840  unsigned Nsample=10;
1841 
1842 
1843  //======================================================================
1850  // If true: populate the bin structure represented by
1860  //======================================================================
1861  void bin_helper(const Vector<Vector<double> >& ray_vertex,
1862  const bool& populate_bin,
1863  Vector<std::pair<unsigned,unsigned> >& intersected_bin,
1864  FiniteElement* el_pt,
1865  std::ofstream& outfile)
1866  {
1867 
1868  // Actually plot
1869  bool plot_it=false;
1870  if (outfile.is_open())
1871  {
1872  if (populate_bin)
1873  {
1874  //Issue a warning
1875  OomphLibWarning(
1876  "Not outputting while bin is being populated...\n",
1877  OOMPH_CURRENT_FUNCTION,
1878  OOMPH_EXCEPTION_LOCATION);
1879  }
1880  else
1881  {
1882  plot_it=true;
1883  }
1884  }
1885 
1886  // Step along ray in increasing y direction
1887  Vector<double> lower_ray_vertex=ray_vertex[0];
1888  Vector<double> upper_ray_vertex=ray_vertex[1];
1889  if (ray_vertex[1][1]<ray_vertex[0][1])
1890  {
1891  lower_ray_vertex=ray_vertex[1];
1892  upper_ray_vertex=ray_vertex[0];
1893  }
1894 
1895  // Output
1896  if (plot_it)
1897  {
1898  outfile
1899  << "ZONE T=\"ray\"\n"
1900  <<ray_vertex[0][0] << " " << ray_vertex[0][1] << "\n"
1901  <<ray_vertex[1][0] << " " << ray_vertex[1][1] << "\n";
1902  }
1903 
1904  // Add bin that contains the start point of the ray
1905  unsigned ix_start=unsigned((lower_ray_vertex[0]-Min_coord[0])/
1906  (Max_coord[0]-Min_coord[0])*
1907  double(Nx_bin));
1908  unsigned iy_start=unsigned((lower_ray_vertex[1]-Min_coord[1])/
1909  (Max_coord[1]-Min_coord[1])*
1910  double(Ny_bin));
1911  if (populate_bin)
1912  {
1913  Element_in_bin[ix_start][iy_start].insert(el_pt);
1914  }
1915  else
1916  {
1917  intersected_bin.push_back(std::make_pair(ix_start,iy_start));
1918  }
1919 
1920  // Vertical ray?
1921  if (lower_ray_vertex[0]==upper_ray_vertex[0])
1922  {
1923  // Add all the bins above the initial one until end point
1924  unsigned iy_end=unsigned((upper_ray_vertex[1]-Min_coord[1])/
1925  (Max_coord[1]-Min_coord[1])*
1926  double(Ny_bin));
1927  for (unsigned i=iy_start+1;i<=iy_end;i++)
1928  {
1929  if (populate_bin)
1930  {
1931  Element_in_bin[ix_start][i].insert(el_pt);
1932  }
1933  else
1934  {
1935  intersected_bin.push_back(std::make_pair(ix_start,i));
1936  }
1937  }
1938  }
1939  // Non-vertical ray
1940  else
1941  {
1942  // Get the (finite) slope
1943  double slope=
1944  (upper_ray_vertex[1]-lower_ray_vertex[1])/
1945  (upper_ray_vertex[0]-lower_ray_vertex[0]);
1946 
1947  // Current bin level
1948  unsigned iy=iy_start;
1949 
1950  // Upper and lower bounds of current bin level
1951  double y_this_bin_min=Min_coord[1]+double(iy )*Dx[1];
1952  double y_this_bin_max=Min_coord[1]+double(iy+1)*Dx[1];
1953 
1954  // Keep going until we've moved through all the
1955  // bins beyond the y-level of the ray's end point
1956  int unit_offset=1;
1957  while (y_this_bin_min<upper_ray_vertex[1])
1958  {
1959  // Work out x-coordinate of intersection of ray with
1960  // other boundary of its current y-bin level
1961  double x_intersect=Max_coord[0];
1962  if (lower_ray_vertex[0]>upper_ray_vertex[0])
1963  {
1964  x_intersect=Min_coord[0];
1965  }
1966  if (slope!=0.0)
1967  {
1968  x_intersect=lower_ray_vertex[0]+
1969  (y_this_bin_max-lower_ray_vertex[1])/slope;
1970  }
1971 
1972  // Cut it off if it goes outside the bin structure
1973  if (x_intersect>Max_coord[0]) x_intersect=Max_coord[0] ;
1974  if (x_intersect<Min_coord[0]) x_intersect=Min_coord[0] ;
1975 
1976  // What's the x-bin index of that point
1977  unsigned ix_intersect=unsigned((x_intersect-Min_coord[0])/
1978  (Max_coord[0]-Min_coord[0])*
1979  double(Nx_bin));
1980 
1981  // limits for x bins:
1982  unsigned i_lo=ix_start;
1983  unsigned i_hi=ix_intersect;
1984 
1985  // Swap if going to the left and apply
1986  // unit offset. Equal to one initially because the first
1987  // bin has already been filled; in subsequent rows of
1988  // bins, we add all of them.
1989  if (i_lo>i_hi)
1990  {
1991  i_lo=ix_intersect;
1992  i_hi=ix_start-unit_offset;
1993  }
1994  else
1995  {
1996  i_lo+=unit_offset;
1997  }
1998  unit_offset=0;
1999 
2000  // ...but don't fall off the end...
2001  if (i_hi==Nx_bin) i_hi-=1;
2002 
2003  // Now add all the bins at this y-level
2004  for (unsigned i=i_lo;i<=i_hi;i++)
2005  {
2006  // .. but don't go beyond the end of the ray
2007  bool add_it=true;
2008  if (slope>0.0)
2009  {
2010  double x_left_end_bin =Min_coord[0]+double(i )*Dx[0];
2011  if (x_left_end_bin>upper_ray_vertex[0])
2012  {
2013  add_it=false;
2014  }
2015  }
2016  else if (slope<0.0)
2017  {
2018  double x_right_end_bin=Min_coord[0]+double(i+1)*Dx[0];
2019  if (x_right_end_bin<upper_ray_vertex[0])
2020  {
2021  add_it=false;
2022  }
2023  }
2024  else
2025  {
2026  double x_left_end_bin =Min_coord[0]+double(i )*Dx[0];
2027  double x_right_end_bin=Min_coord[0]+double(i+1)*Dx[0];
2028  if (x_left_end_bin>std::max(upper_ray_vertex[0],lower_ray_vertex[0]))
2029  {
2030  add_it=false;
2031  }
2032  else if (x_right_end_bin<std::min(upper_ray_vertex[0],
2033  lower_ray_vertex[0]))
2034  {
2035  add_it=false;
2036  }
2037  }
2038 
2039  if (add_it)
2040  {
2041  if (populate_bin)
2042  {
2043  Element_in_bin[i][iy].insert(el_pt);
2044  }
2045  else
2046  {
2047  intersected_bin.push_back(std::make_pair(i,iy));
2048  }
2049  }
2050  }
2051 
2052  // Now update the counters: x bin level starts at the
2053  // point of intersection with the next level...
2054  ix_start=ix_intersect;
2055 
2056  // ...upper and lower boundary moves one level up
2057  y_this_bin_min+=Dx[1];
2058  y_this_bin_max+=Dx[1];
2059 
2060  // ...and bump up the level itself
2061  iy++;
2062  }
2063  }
2064 
2065  // Plot the intersected bins?
2066  if (plot_it)
2067  {
2068  unsigned nbin=intersected_bin.size();
2069  for (unsigned i=0;i<nbin;i++)
2070  {
2071  unsigned ix=intersected_bin[i].first;
2072  unsigned iy=intersected_bin[i].second;
2073 
2074  double x_lo=Min_coord[0]+double(ix)*Dx[0];
2075  double x_hi=x_lo+Dx[0];
2076  double y_lo=Min_coord[1]+double(iy)*Dx[1];
2077  double y_hi=y_lo+Dx[1];
2078 
2079  outfile << "ZONE\n"
2080  << x_lo << " " << y_lo << "\n"
2081  << x_hi << " " << y_lo << "\n"
2082  << x_hi << " " << y_hi << "\n"
2083  << x_lo << " " << y_hi << "\n"
2084  << x_lo << " " << y_lo << "\n";
2085  }
2086  }
2087  }
2088 
2089 
2090  //=================================================================
2091  // Doc populated bins
2092  //=================================================================
2093  void doc_bins(std::ofstream& bin_file)
2094  {
2095  // Loop over bins
2096  for (unsigned ix=0;ix<Nx_bin;ix++)
2097  {
2098  for (unsigned iy=0;iy<Ny_bin;iy++)
2099  {
2100  // Anybody at home?
2101  if (Element_in_bin[ix][iy].size()!=0)
2102  {
2103  double x_lo=Min_coord[0]+double(ix)*Dx[0];
2104  double x_hi=x_lo+Dx[0];
2105  double y_lo=Min_coord[1]+double(iy)*Dx[1];
2106  double y_hi=y_lo+Dx[1];
2107 
2108  bin_file << "ZONE I=2, J=2\n"
2109  << x_lo << " " << y_lo << "\n"
2110  << x_hi << " " << y_lo << "\n"
2111  << x_lo << " " << y_hi << "\n"
2112  << x_hi << " " << y_hi << "\n";
2113 
2114  }
2115  }
2116  }
2117  }
2118 
2119 
2120 
2121  //=================================================================
2122  // Doc sample points of Stefan Boltzmann elements
2123  //=================================================================
2124  void doc_sample_points(std::ofstream& outfile,
2125  const Vector<FiniteElement*>& sb_face_element_pt)
2126  {
2127  // Vector for coordinates of sample point
2128  Vector<double> sample_point(2);
2129  Vector<double> s(1);
2130 
2131  // Loop over all face elements
2132  unsigned nel=sb_face_element_pt.size();
2133  for (unsigned e=0;e<nel;e++)
2134  {
2135  StefanBoltzmannRadiationBase* el_pt=
2136  dynamic_cast<StefanBoltzmannRadiationBase*>(
2137  sb_face_element_pt[e]);
2138 #ifdef PARANOID
2139  if (el_pt==0)
2140  {
2141  std::stringstream error_message;
2142  error_message << "Failed to cast possible intersecting element "
2143  << e
2144  << " to StefanBoltzmannRadiationBase";
2145  throw OomphLibError(
2146  error_message.str(),
2147  OOMPH_CURRENT_FUNCTION,
2148  OOMPH_EXCEPTION_LOCATION);
2149  }
2150 #endif
2151 
2152  // Loop over sampling points
2153  for (unsigned j_sample=0;j_sample<Nsample;j_sample++)
2154  {
2155  // Get vector of local coordinates of plot point j
2156  // as second vertex
2157  el_pt->get_s_plot(j_sample,Nsample,s);
2158 
2159  // Get coordinate of vertex
2160  el_pt->interpolated_x(s,sample_point);
2161 
2162  outfile << sample_point[0] << " " << sample_point[1] << "\n";
2163  }
2164  }
2165  }
2166 
2167 
2168  //=======================================================================
2172  //=======================================================================
2173  void setup_stefan_boltzmann_visibility(const Vector<FiniteElement*>&
2174  sb_face_element_pt)
2175  {
2176  // Output file for debugging
2177  const bool plot_it=false;
2178  std::ofstream some_file;
2179  char filename[100];
2180 
2181  // Loop over all face elements to wipe previous information
2182  unsigned nel=sb_face_element_pt.size();
2183  for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
2184  {
2185  StefanBoltzmannRadiationBase* illuminated_el_pt=
2186  dynamic_cast<StefanBoltzmannRadiationBase*>(
2187  sb_face_element_pt[e_illuminated]);
2188 #ifdef PARANOID
2189  if (illuminated_el_pt==0)
2190  {
2191  std::stringstream error_message;
2192  error_message << "Failed to cast illuminated element "
2193  << e_illuminated
2194  << " to StefanBoltzmannRadiationBase";
2195  throw OomphLibError(
2196  error_message.str(),
2197  OOMPH_CURRENT_FUNCTION,
2198  OOMPH_EXCEPTION_LOCATION);
2199  }
2200 #endif
2201 
2202  // Forget any previous information
2203  illuminated_el_pt->wipe_stefan_boltzmann_illumination_info();
2204  }
2205 
2206 
2207  double t_start=TimingHelpers::timer();
2208 
2209  // Vector to illuminated/ing Gauss points
2210  Vector<double> r_illuminated(2);
2211  Vector<double> s_illuminated(1);
2212  Vector<double> unit_normal_illuminated(2);
2213  Vector<double> r_illuminating(2);
2214  Vector<double> s_illuminating(1);
2215  Vector<double> unit_normal_illuminating(2);
2216 
2217  // Setup bin structure
2218  const bool use_bins=true;
2219  if (use_bins)
2220  {
2221  double t_start_bin=TimingHelpers::timer();
2222 
2223  // Each bin stores FiniteElements that have any sample point in it
2224  Element_in_bin.clear();
2225  Element_in_bin.resize(Nx_bin);
2226  for (unsigned i=0;i<Nx_bin;i++)
2227  {
2228  Element_in_bin[i].resize(Ny_bin);
2229  }
2230 
2231 
2232  // Reset max/min coords
2233  Max_coord.clear();
2234  Max_coord.resize(2,-DBL_MAX);
2235  Min_coord.clear();
2236  Min_coord.resize(2, DBL_MAX);
2237 
2238  // Initial sweep: Get max/min coords
2239  //----------------------------------
2240 
2241  // Loop over all face elements
2242  unsigned nel=sb_face_element_pt.size();
2243  for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
2244  {
2245  StefanBoltzmannRadiationBase* illuminated_el_pt=
2246  dynamic_cast<StefanBoltzmannRadiationBase*>(
2247  sb_face_element_pt[e_illuminated]);
2248 #ifdef PARANOID
2249  if (illuminated_el_pt==0)
2250  {
2251  std::stringstream error_message;
2252  error_message << "Failed to cast illuminated element "
2253  << e_illuminated
2254  << " to StefanBoltzmannRadiationBase";
2255  throw OomphLibError(
2256  error_message.str(),
2257  OOMPH_CURRENT_FUNCTION,
2258  OOMPH_EXCEPTION_LOCATION);
2259  }
2260 #endif
2261 
2262  // Loop over iluminated integration points
2263  unsigned nint=illuminated_el_pt->integral_pt()->nweight();
2264  for (unsigned ipt_illuminated=0;ipt_illuminated<nint;ipt_illuminated++)
2265  {
2266  // Local coordinate of integration point
2267  s_illuminated[0]=
2268  illuminated_el_pt->integral_pt()->knot(ipt_illuminated,0);
2269 
2270  // No recycling of shape fcts -- may as well call interpolated_x
2271  // directly
2272  illuminated_el_pt->interpolated_x(s_illuminated,r_illuminated);
2273 
2274  for (unsigned i=0;i<2;i++)
2275  {
2276  if (r_illuminated[i]>Max_coord[i]) Max_coord[i]=r_illuminated[i];
2277  if (r_illuminated[i]<Min_coord[i]) Min_coord[i]=r_illuminated[i];
2278  }
2279  }
2280  }
2281 
2282  // Allow for offset
2283  double percentage_offset=5.0;
2284  for (unsigned i=0;i<2;i++)
2285  {
2286  double dx=Max_coord[i]-Min_coord[i];
2287  Max_coord[i]=Max_coord[i]+percentage_offset/100.0*dx;
2288  Min_coord[i]=Min_coord[i]-percentage_offset/100.0*dx;
2289  }
2290 
2291  // Update increments
2292  Dx[0]=(Max_coord[0]-Min_coord[0])/double(Nx_bin);
2293  Dx[1]=(Max_coord[1]-Min_coord[1])/double(Ny_bin);
2294 
2295 
2296  const bool test_horizontal_and_vertical_rays=false;
2297  if (test_horizontal_and_vertical_rays)
2298  {
2299 
2300  // Test for left to right horizontal ray
2301  {
2302  // Populate bin by associating all bins intersected by
2303  // ray from one to the other sample vertex with the
2304  // element (note that the ray can span multiple bins!)
2305  Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
2306  bool populate_bin=true;
2307  std::ofstream dummy_file;
2308  Vector<Vector<double> > sample_vertex(2);
2309  sample_vertex[0].resize(2);
2310  sample_vertex[1].resize(2);
2311  sample_vertex[0][0]=-0.9;
2312  sample_vertex[0][1]=0.9;
2313  sample_vertex[1][0]=0.9;
2314  sample_vertex[1][1]=0.9;
2315 
2316  StefanBoltzmannRadiationBase* el_pt=
2317  dynamic_cast<StefanBoltzmannRadiationBase*>(
2318  sb_face_element_pt[0]);
2319 
2320  bin_helper(sample_vertex,
2321  populate_bin,
2322  dummy_intersected_bin,
2323  el_pt,
2324  dummy_file);
2325  }
2326 
2327 
2328 
2329 
2330  // Test for right to left horizontal ray
2331  {
2332  // Populate bin by associating all bins intersected by
2333  // ray from one to the other sample vertex with the
2334  // element (note that the ray can span multiple bins!)
2335  Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
2336  bool populate_bin=true;
2337  std::ofstream dummy_file;
2338  Vector<Vector<double> > sample_vertex(2);
2339  sample_vertex[0].resize(2);
2340  sample_vertex[1].resize(2);
2341  sample_vertex[0][0]= 0.9;
2342  sample_vertex[0][1]=-0.9;
2343  sample_vertex[1][0]=-0.9;
2344  sample_vertex[1][1]=-0.9;
2345 
2346  StefanBoltzmannRadiationBase* el_pt=
2347  dynamic_cast<StefanBoltzmannRadiationBase*>(
2348  sb_face_element_pt[0]);
2349 
2350  bin_helper(sample_vertex,
2351  populate_bin,
2352  dummy_intersected_bin,
2353  el_pt,
2354  dummy_file);
2355  }
2356 
2357 
2358  // Test for bottom to top vertical ray
2359  {
2360  // Populate bin by associating all bins intersected by
2361  // ray from one to the other sample vertex with the
2362  // element (note that the ray can span multiple bins!)
2363  Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
2364  bool populate_bin=true;
2365  std::ofstream dummy_file;
2366  Vector<Vector<double> > sample_vertex(2);
2367  sample_vertex[0].resize(2);
2368  sample_vertex[1].resize(2);
2369  sample_vertex[0][0]=-0.7;
2370  sample_vertex[0][1]=-0.9;
2371  sample_vertex[1][0]=-0.7;
2372  sample_vertex[1][1]= 0.9;
2373 
2374  StefanBoltzmannRadiationBase* el_pt=
2375  dynamic_cast<StefanBoltzmannRadiationBase*>(
2376  sb_face_element_pt[0]);
2377 
2378  bin_helper(sample_vertex,
2379  populate_bin,
2380  dummy_intersected_bin,
2381  el_pt,
2382  dummy_file);
2383  }
2384 
2385 
2386  // Test for to to bottom vertical ray
2387  {
2388  // Populate bin by associating all bins intersected by
2389  // ray from one to the other sample vertex with the
2390  // element (note that the ray can span multiple bins!)
2391  Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
2392  bool populate_bin=true;
2393  std::ofstream dummy_file;
2394  Vector<Vector<double> > sample_vertex(2);
2395  sample_vertex[0].resize(2);
2396  sample_vertex[1].resize(2);
2397  sample_vertex[0][0]= 0.7;
2398  sample_vertex[0][1]= 0.9;
2399  sample_vertex[1][0]= 0.7;
2400  sample_vertex[1][1]=-0.9;
2401 
2402  StefanBoltzmannRadiationBase* el_pt=
2403  dynamic_cast<StefanBoltzmannRadiationBase*>(
2404  sb_face_element_pt[0]);
2405 
2406  bin_helper(sample_vertex,
2407  populate_bin,
2408  dummy_intersected_bin,
2409  el_pt,
2410  dummy_file);
2411  }
2412 
2413  } // end of test for horizontal and vertical rays
2414 
2415  // Setup bin for checking intersections with rays
2416  //-----------------------------------------------
2417 
2418  // Loop over all face elements
2419  nel=sb_face_element_pt.size();
2420  for (unsigned e=0;e<nel;e++)
2421  {
2422  StefanBoltzmannRadiationBase* el_pt=
2423  dynamic_cast<StefanBoltzmannRadiationBase*>(
2424  sb_face_element_pt[e]);
2425 #ifdef PARANOID
2426  if (el_pt==0)
2427  {
2428  std::stringstream error_message;
2429  error_message << "Failed to cast possible intersecting element "
2430  << e
2431  << " to StefanBoltzmannRadiationBase";
2432  throw OomphLibError(
2433  error_message.str(),
2434  OOMPH_CURRENT_FUNCTION,
2435  OOMPH_EXCEPTION_LOCATION);
2436  }
2437 #endif
2438 
2439  // Loop over ALL sampling points along element
2440  Vector<Vector<double> > sample_vertex(2);
2441  sample_vertex[0].resize(2);
2442  sample_vertex[1].resize(2);
2443 
2444  // Get vector of local coordinates of plot point j
2445  // as first vertex
2446  Vector<double> s(1);
2447  unsigned j_sample=0;
2448  el_pt->get_s_plot(j_sample,Nsample,s);
2449 
2450  // Get coordinate of first vertex
2451  el_pt->interpolated_x(s,sample_vertex[0]);
2452 
2453  // Loop over remaining sampling points
2454  for (unsigned j_sample=1;j_sample<Nsample;j_sample++)
2455  {
2456  // Get vector of local coordinates of plot point j
2457  // as second vertex
2458  el_pt->get_s_plot(j_sample,Nsample,s);
2459 
2460  // Get coordinate of vertex
2461  el_pt->interpolated_x(s,sample_vertex[1]);
2462 
2463  // Populate bin by associating all bins intersected by
2464  // ray from one to the other sample vertex with the
2465  // element (note that the ray can span multiple bins!)
2466  Vector<std::pair<unsigned,unsigned> > dummy_intersected_bin;
2467  bool populate_bin=true;
2468  std::ofstream dummy_file;
2469  bin_helper(sample_vertex,
2470  populate_bin,
2471  dummy_intersected_bin,
2472  el_pt,
2473  dummy_file);
2474 
2475  // Shift along
2476  sample_vertex[0]=sample_vertex[1];
2477  }
2478  }
2479 
2480  double t_end_bin=TimingHelpers::timer();
2481  oomph_info << "Time for setting up bin: "
2482  << t_end_bin-t_start_bin << std::endl;
2483 
2484  // Number of elements
2485  nel=sb_face_element_pt.size();
2486 
2487  // Assume all elements have the same number of integration points
2488  unsigned nintpt_all=sb_face_element_pt[0]->integral_pt()->nweight();
2489 
2490  // Exploit symmetry during setup
2491  const bool exploit_symmetry=true;
2492 
2493  // Helper lookup scheme to exploit symmetry of interaction:
2494  // If integration point i_ed in element e_ed is illuminatED by
2495  // integration point i_ing in element e_ing (which is the
2496  // illuminatING one!) then the reverse is true too
2497  Vector<Vector<Vector<Vector<unsigned> > > > aux;
2498  if (exploit_symmetry)
2499  {
2500  nintpt_all=sb_face_element_pt[0]->integral_pt()->nweight();
2501  aux.resize(nel);
2502  for (unsigned e=0;e<nel;e++)
2503  {
2504  aux[e].resize(nintpt_all);
2505  for (unsigned ipt=0;ipt<nintpt_all;ipt++)
2506  {
2507  aux[e][ipt].resize(nel);
2508  }
2509  }
2510  }
2511 
2512  // Loop over all face elements -- viewed as illuminated ones
2513  for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
2514  {
2515  StefanBoltzmannRadiationBase* illuminated_el_pt=
2516  dynamic_cast<StefanBoltzmannRadiationBase*>(
2517  sb_face_element_pt[e_illuminated]);
2518 
2519 #ifdef PARANOID
2520  if (illuminated_el_pt==0)
2521  {
2522  std::stringstream error_message;
2523  error_message << "Failed to cast illuminated element "
2524  << e_illuminated
2525  << " to StefanBoltzmannRadiationBase";
2526  throw OomphLibError(
2527  error_message.str(),
2528  OOMPH_CURRENT_FUNCTION,
2529  OOMPH_EXCEPTION_LOCATION);
2530  }
2531 #endif
2532 
2533  // Recast to FaceElement
2534  FaceElement* illuminated_face_el_pt=
2535  dynamic_cast<FaceElement*>(illuminated_el_pt);
2536 
2537  // Loop over iluminated integration points
2538  unsigned nint=illuminated_el_pt->integral_pt()->nweight();
2539  for (unsigned ipt_illuminated=0;ipt_illuminated<nint;ipt_illuminated++)
2540  {
2541  // Local coordinate of integration point
2542  s_illuminated[0]=
2543  illuminated_el_pt->integral_pt()->knot(ipt_illuminated,0);
2544 
2545  // No recycling of shape fcts -- may as well call interpolated_x
2546  // directly
2547  illuminated_el_pt->interpolated_x(s_illuminated,r_illuminated);
2548 
2549  // Get outer unit normal
2550  illuminated_face_el_pt->outer_unit_normal(s_illuminated,
2551  unit_normal_illuminated);
2552 
2553  // Now loop over all other elements (the iluminating ones)
2554  // Note: this includes the current one because its integration
2555  // points may illuminate each other!
2556  unsigned e_lo=0;
2557  if (exploit_symmetry) e_lo=e_illuminated;
2558  for (unsigned e_illuminating=e_lo;e_illuminating<nel;e_illuminating++)
2559  {
2560  StefanBoltzmannRadiationBase* illuminating_el_pt=
2561  dynamic_cast<StefanBoltzmannRadiationBase*>(
2562  sb_face_element_pt[e_illuminating]);
2563 
2564 #ifdef PARANOID
2565  if (illuminating_el_pt==0)
2566  {
2567  std::stringstream error_message;
2568  error_message << "Failed to cast illuminating element "
2569  << e_illuminating
2570  << " to StefanBoltzmannRadiationBase";
2571  throw OomphLibError(
2572  error_message.str(),
2573  OOMPH_CURRENT_FUNCTION,
2574  OOMPH_EXCEPTION_LOCATION);
2575  }
2576 #endif
2577 
2578  // Recast to FaceElement
2579  FaceElement* illuminating_face_el_pt=
2580  dynamic_cast<FaceElement*>(illuminating_el_pt);
2581 
2582  // Storage to accumulate indices of integration points in
2583  // illuminating face element that is visible from
2584  // current integration point
2585  Vector<unsigned> illuminating_integration_point_index;
2586 
2587  // Loop over iluminating integration points
2588  unsigned nint=illuminating_el_pt->integral_pt()->nweight();
2589  for (unsigned ipt_illuminating=0;ipt_illuminating<nint;
2590  ipt_illuminating++)
2591  {
2592  // Local coordinate of integration point
2593  s_illuminating[0]=
2594  illuminating_el_pt->integral_pt()->knot(ipt_illuminating,0);
2595 
2596  // No recycling of shape fcts -- may as well call interpolated_x
2597  // directly
2598  illuminating_el_pt->interpolated_x(s_illuminating,r_illuminating);
2599 
2600  // Get outer unit normal
2601  illuminating_face_el_pt->
2602  outer_unit_normal(s_illuminating,
2603  unit_normal_illuminating);
2604 
2605  // Get vector from illuminating point to illuminated point
2606  Vector<double> ray(2);
2607  ray[0]=r_illuminated[0]-r_illuminating[0];
2608  ray[1]=r_illuminated[1]-r_illuminating[1];
2609 
2610 
2611  Vector<Vector<double> > ray_vertex(2);
2612  ray_vertex[0].resize(2);
2613  ray_vertex[1].resize(2);
2614  ray_vertex[0][0]=r_illuminated[0];
2615  ray_vertex[0][1]=r_illuminated[1];
2616  ray_vertex[1][0]=r_illuminating[0];
2617  ray_vertex[1][1]=r_illuminating[1];
2618 
2619  // Do we radiate away from the positive face of the
2620  // illuminating point?
2621  double dot_illuminating=
2622  (unit_normal_illuminating[0]*ray[0]+
2623  unit_normal_illuminating[1]*ray[1]);
2624  if (dot_illuminating>0.0)
2625  {
2626 
2627  // Do we radiate onto from the positive face of the
2628  // illuminated point?
2629  double dot_illuminated=
2630  (unit_normal_illuminated[0]*ray[0]+
2631  unit_normal_illuminated[1]*ray[1]);
2632  if (dot_illuminated<0.0)
2633  {
2634 
2635  // Does the ray (a finite-length segment) from
2636  // radiating to radiated point intersect any other elements?
2637 
2638  // Vector of pairs of bin coordinates that
2639  // intersect ray
2640  Vector<std::pair<unsigned,unsigned> > intersected_bin;
2641  if (use_bins)
2642  {
2643  if (plot_it)
2644  {
2645  sprintf(filename,"RESLT/latest_ray.dat");
2646  some_file.open(filename);
2647  }
2648 
2649  // Find bins that are intersected by ray
2650  Vector<std::pair<unsigned,unsigned> > intersected_bin;
2651  bool populate_bin=false;
2652  FiniteElement* dummy_el_pt=0;
2653  bin_helper(ray_vertex,
2654  populate_bin,
2655  intersected_bin,
2656  dummy_el_pt,
2657  some_file);
2658 
2659  // End plot
2660  if (plot_it)
2661  {
2662  some_file.close();
2663  pause("done latest ray");
2664  }
2665 
2666  // Storage for vertices of possible intersection with ray
2667  Vector<Vector<double> > segment_vertex(2);
2668  segment_vertex[0].resize(2);
2669  segment_vertex[1].resize(2);
2670 
2671  // Search through all elements in bins along ray
2672  bool have_intersection=false;
2673  unsigned nbin=intersected_bin.size();
2674  for (unsigned b=0;b<nbin;b++)
2675  {
2676  // Get bin indices
2677  unsigned i=intersected_bin[b].first;
2678  unsigned j=intersected_bin[b].second;
2679 
2680  // Loop over elements in that bin
2681  for (std::set<FiniteElement*>::iterator it=
2682  Element_in_bin[i][j].begin();
2683  it!=Element_in_bin[i][j].end();it++)
2684  {
2685  // Get possibly blocking element
2686  StefanBoltzmannRadiationBase* el_pt=
2687  dynamic_cast<StefanBoltzmannRadiationBase*>(*it);
2688 
2689 #ifdef PARANOID
2690  if (illuminated_el_pt==0)
2691  {
2692  std::stringstream error_message;
2693  error_message
2694  << "Failed to cast possibly blocking element "
2695  << " to StefanBoltzmannRadiationBase";
2696  throw OomphLibError(
2697  error_message.str(),
2698  OOMPH_CURRENT_FUNCTION,
2699  OOMPH_EXCEPTION_LOCATION);
2700  }
2701 #endif
2702  // Skip intersections with illuming/illuminated element
2703  if ((el_pt!=illuminated_el_pt)&&
2704  (el_pt!=illuminating_el_pt))
2705  {
2706  // Loop over sampling points along element, only
2707  // up to Nsample-1 because we're forming segment
2708  // with current and next sampling point.
2709  for (unsigned j_sample=0;j_sample<Nsample-1;j_sample++)
2710  {
2711  // Get cector of local coordinates of plot point j
2712  // as first vertex
2713  Vector<double> s(1);
2714  el_pt->get_s_plot(j_sample,Nsample,s);
2715 
2716  // Get coordinate of first vertex
2717  el_pt->interpolated_x(s,segment_vertex[0]);
2718 
2719  // Get cector of local coordinates of plot point j+1
2720  // as second vertex
2721  el_pt->get_s_plot(j_sample+1,Nsample,s);
2722 
2723  // Get coordinate of second vertex
2724  el_pt->interpolated_x(s,segment_vertex[1]);
2725 
2726  // Check intersection
2727  bool intersection=IntersectionChecker::intersects(
2728  segment_vertex,ray_vertex);
2729 
2730  // Bail out
2731  if (intersection)
2732  {
2733  have_intersection=true;
2734  break;
2735  }
2736  } // end of loop over sampling points
2737  } // endif for self-intersection
2738  if (have_intersection) break;
2739  } // end of loop over elements in bin
2740  if (have_intersection) break;
2741  }// End of loop over bins that contain ray
2742 
2743 
2744  // No intersection
2745  if (!have_intersection)
2746  {
2747  // Visible, so add integration point
2748  illuminating_integration_point_index.
2749  push_back(ipt_illuminating);
2750  }
2751  }
2752  // No bins: Brute force search loop
2753  else
2754  {
2755  // Intersection for star-shaped (non-convex) polygon
2756  // can only be detected in O(n^2) operations.
2757  Vector<Vector<double> > segment_vertex(2);
2758  segment_vertex[0].resize(2);
2759  segment_vertex[1].resize(2);
2760 
2761  // Loop over all segments to test for intersection
2762  bool have_intersection=false;
2763  for (unsigned e_intersect=0;e_intersect<nel;
2764  e_intersect++)
2765  {
2766 
2767  // Skip intersection with illuminating/ed elements
2768  if (! ( (e_intersect==e_illuminated) ||
2769  (e_intersect==e_illuminating) ) )
2770  {
2771 
2772  // Loop over sampling points along element
2773  for (unsigned j=0;j<Nsample-1;j++)
2774  {
2775  // Get cector of local coordinates of plot point j
2776  // as first vertex
2777  Vector<double> s(1);
2778  sb_face_element_pt[e_intersect]->
2779  get_s_plot(j,Nsample,s);
2780 
2781  // Get coordinate of first vertex
2782  sb_face_element_pt[e_intersect]->
2783  interpolated_x(s,segment_vertex[0]);
2784 
2785  // Get cector of local coordinates of plot point j+1
2786  // as second vertex
2787  sb_face_element_pt[e_intersect]->
2788  get_s_plot(j+1,Nsample,s);
2789 
2790  // Get coordinate of second vertex
2791  sb_face_element_pt[e_intersect]->
2792  interpolated_x(s,segment_vertex[1]);
2793 
2794  // Check intersection
2795  bool intersection=IntersectionChecker::intersects(
2796  segment_vertex,ray_vertex);
2797 
2798  if (intersection)
2799  {
2800  have_intersection=true;
2801  break;
2802  }
2803  }
2804  }
2805  if (have_intersection)
2806  {
2807  break;
2808  }
2809  }
2810 
2811  // No intersection
2812  if (!have_intersection)
2813  {
2814  // Visible, so add integration point
2815  illuminating_integration_point_index.
2816  push_back(ipt_illuminating);
2817  }
2818 
2819  } // end if for brute force (rather than bin-based) search loop
2820 
2821  }
2822  }
2823  }
2824 
2825  // Done all possibly illuminating integration points in
2826  // current possibly illuminating element
2827  unsigned npt=illuminating_integration_point_index.size();
2828  if (npt>0)
2829  {
2830  // Add info
2831  illuminated_el_pt->add_stefan_boltzmann_illumination_info(
2832  ipt_illuminated,illuminating_el_pt,
2833  illuminating_integration_point_index);
2834 
2835  // Set up reverse scheme
2836  if (exploit_symmetry)
2837  {
2838  for (unsigned i_ing=0;i_ing<npt;i_ing++)
2839  {
2840  aux[e_illuminating]
2841  [illuminating_integration_point_index[i_ing]]
2842  [e_illuminated].push_back(ipt_illuminated);
2843  }
2844  }
2845  }
2846  }
2847  }
2848  }
2849 
2850 
2851  // Now do other half of dependencies
2852  //----------------------------------
2853  if (exploit_symmetry)
2854  {
2855  // Loop over all face elements -- viewed as illuminated ones
2856  for (unsigned e_illuminated=0;e_illuminated<nel;e_illuminated++)
2857  {
2858  StefanBoltzmannRadiationBase* illuminated_el_pt=
2859  dynamic_cast<StefanBoltzmannRadiationBase*>(
2860  sb_face_element_pt[e_illuminated]);
2861 
2862 #ifdef PARANOID
2863  if (illuminated_el_pt==0)
2864  {
2865  std::stringstream error_message;
2866  error_message << "Failed to cast illuminated element "
2867  << e_illuminated
2868  << " to StefanBoltzmannRadiationBase";
2869  throw OomphLibError(
2870  error_message.str(),
2871  OOMPH_CURRENT_FUNCTION,
2872  OOMPH_EXCEPTION_LOCATION);
2873  }
2874 #endif
2875 
2876  // Loop over iluminated integration points
2877  unsigned nint=illuminated_el_pt->integral_pt()->nweight();
2878  for (unsigned ipt_illuminated=0;ipt_illuminated<nint;ipt_illuminated++)
2879  {
2880  // Now loop over all other elements (the iluminating ones)
2881  for (unsigned e_illuminating=0;
2882  e_illuminating<e_illuminated;e_illuminating++)
2883  {
2884  StefanBoltzmannRadiationBase* illuminating_el_pt=
2885  dynamic_cast<StefanBoltzmannRadiationBase*>(
2886  sb_face_element_pt[e_illuminating]);
2887 
2888 #ifdef PARANOID
2889  if (illuminating_el_pt==0)
2890  {
2891  std::stringstream error_message;
2892  error_message << "Failed to cast illuminating element "
2893  << e_illuminating
2894  << " to StefanBoltzmannRadiationBase";
2895  throw OomphLibError(
2896  error_message.str(),
2897  OOMPH_CURRENT_FUNCTION,
2898  OOMPH_EXCEPTION_LOCATION);
2899  }
2900 #endif
2901  // Get illuminating integration points
2902  Vector<unsigned> illuminating_integration_point_index=
2903  aux[e_illuminated][ipt_illuminated][e_illuminating];
2904  unsigned npt=illuminating_integration_point_index.size();
2905  if (npt>0)
2906  {
2907  // Add info
2908  illuminated_el_pt->add_stefan_boltzmann_illumination_info(
2909  ipt_illuminated,illuminating_el_pt,
2910  illuminating_integration_point_index);
2911  }
2912  }
2913  }
2914  }
2915  }
2916 
2917  double t_end=TimingHelpers::timer();
2918  oomph_info << "Time for setting up mutual Stefan Boltzmann radiation: "
2919  << t_end-t_start << std::endl;
2920  }
2921 
2922  }
2923 
2924 }
2925 }
2926 
2927 
2928 
2929 #endif
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
atan2
AnnoyingScalar atan2(const AnnoyingScalar &y, const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:139
sin
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::DShape
Definition: shape.h:278
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
oomph::FaceElement
Definition: elements.h:4338
oomph::FaceElement::face_index
int & face_index()
Definition: elements.h:4626
oomph::FaceElement::outer_unit_normal
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
oomph::FaceElement::zeta_nodal
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Definition: elements.h:4497
oomph::FaceElement::interpolated_x
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Definition: elements.h:4528
oomph::FaceGeometry
Definition: elements.h:4998
oomph::FiniteElement
Definition: elements.h:1313
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::size
double size() const
Definition: elements.cc:4290
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::FiniteElement::dshape_local_at_knot
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Definition: elements.cc:3239
oomph::FiniteElement::dim
unsigned dim() const
Definition: elements.h:2611
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::get_s_plot
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Definition: elements.h:3148
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Definition: elements.h:2317
oomph::FiniteElement::nodal_dimension
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
oomph::FiniteElement::set_integration_scheme
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3210
oomph::GeneralisedElement::add_external_data
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:307
oomph::Integral
Definition: integral.h:49
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::Node
Definition: nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
oomph::Node::value
double value(const unsigned &i) const
Definition: nodes.cc:2408
oomph::OomphLibError
Definition: oomph_definitions.h:222
oomph::OomphLibWarning
Definition: oomph_definitions.h:267
oomph::RefineableElement
Definition: refineable_elements.h:97
oomph::Shape
Definition: shape.h:76
oomph::SolarRadiationBase
Definition: temporary_stefan_boltzmann_elements.h:150
oomph::SolarRadiationBase::output_limiting_angles
void output_limiting_angles(std::ostream &outfile)
Output illumination angles for all integration points.
Definition: temporary_stefan_boltzmann_elements.h:1020
oomph::SolarRadiationBase::smoothed_sun_shadow_is_enabled
bool smoothed_sun_shadow_is_enabled()
Disable smoothing of shadow.
Definition: temporary_stefan_boltzmann_elements.h:194
oomph::SolarRadiationBase::Diffuse_limit_angles
Vector< std::pair< double, double > > Diffuse_limit_angles
Definition: temporary_stefan_boltzmann_elements.h:270
oomph::SolarRadiationBase::output_atmospheric_radiation
void output_atmospheric_radiation(std::ostream &outfile)
Output diffuse and direct radiation.
Definition: temporary_stefan_boltzmann_elements.h:1068
oomph::SolarRadiationBase::output_diffuse_radiation_cone_max_angle
void output_diffuse_radiation_cone_max_angle(std::ostream &outfile, const double &radius)
Definition: temporary_stefan_boltzmann_elements.h:920
oomph::SolarRadiationBase::disable_smoothed_sun_shadow
void disable_smoothed_sun_shadow()
Disable smoothing of shadow.
Definition: temporary_stefan_boltzmann_elements.h:188
oomph::SolarRadiationBase::atmospheric_radiation_fct_pt
void(*&)(const double &time, double &solar_flux_magnitude, Vector< double > &solar_flux_unit_vector, double &total_diffuse_radiation) atmospheric_radiation_fct_pt()
Reference to the atmospheric radiation function pointer.
Definition: temporary_stefan_boltzmann_elements.h:209
oomph::SolarRadiationBase::enable_smoothed_sun_shadow
void enable_smoothed_sun_shadow(const double &alpha_tanh_smooth_sun_shadow=100.0)
Definition: temporary_stefan_boltzmann_elements.h:180
oomph::SolarRadiationBase::SolarRadiationBase
SolarRadiationBase()
Constructor.
Definition: temporary_stefan_boltzmann_elements.h:155
oomph::SolarRadiationBase::output_diffuse_radiation_cone
void output_diffuse_radiation_cone(std::ostream &outfile, const double &radius)
Output cone of diffuse radiation for all integration points.
Definition: temporary_stefan_boltzmann_elements.h:865
oomph::SolarRadiationBase::Smoothed_sun_shadow
bool Smoothed_sun_shadow
Use tanh profile to smooth solar shadows.
Definition: temporary_stefan_boltzmann_elements.h:261
oomph::SolarRadiationBase::update_limiting_angles
void update_limiting_angles(const Vector< Node * > &shielding_node_pt)
Definition: temporary_stefan_boltzmann_elements.h:528
oomph::SolarRadiationBase::atmospheric_radiation
virtual double atmospheric_radiation(const unsigned &intpt, const double &time, const Vector< double > &x, const Vector< double > &n)
Definition: temporary_stefan_boltzmann_elements.h:309
oomph::SolarRadiationBase::Alpha_tanh_smooth_sun_shadow
double Alpha_tanh_smooth_sun_shadow
Factor for tanh profile to smooth solar shadows.
Definition: temporary_stefan_boltzmann_elements.h:264
oomph::SolarRadiationBase::alpha_tanh_smooth_sun_shadow
double alpha_tanh_smooth_sun_shadow()
Definition: temporary_stefan_boltzmann_elements.h:202
oomph::SolarRadiationBase::Atmospheric_radiation_fct_pt
void(* Atmospheric_radiation_fct_pt)(const double &time, double &solar_flux_magnitude, Vector< double > &solar_flux_unit_vector, double &total_diffuse_radiation)
Definition: temporary_stefan_boltzmann_elements.h:276
oomph::SolarRadiationBase::check_quadrant_jump
void check_quadrant_jump(const double &x_prev, const double &y_prev, const double &x_next, const double &y_next, bool &crossing_quadrants, int &n_winding_increment, std::string &error_string, bool &no_problem)
Definition: temporary_stefan_boltzmann_elements.h:409
oomph::SolarRadiationBase::output_diffuse_radiation_cone_min_angle
void output_diffuse_radiation_cone_min_angle(std::ostream &outfile, const double &radius)
Definition: temporary_stefan_boltzmann_elements.h:969
oomph::SolidNode
Definition: nodes.h:1686
oomph::SolidNode::variable_position_pt
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1765
oomph::StefanBoltzmannRadiationBase
Definition: temporary_stefan_boltzmann_elements.h:1132
oomph::StefanBoltzmannRadiationBase::Stefan_boltzmann_illumination_info
Vector< Vector< std::pair< StefanBoltzmannRadiationBase *, Vector< unsigned > > > > Stefan_boltzmann_illumination_info
Definition: temporary_stefan_boltzmann_elements.h:1361
oomph::StefanBoltzmannRadiationBase::add_stefan_boltzmann_illumination_info
void add_stefan_boltzmann_illumination_info(const unsigned &ipt, StefanBoltzmannRadiationBase *illuminating_el_pt, Vector< unsigned > &illuminating_integration_point_index, const bool &add_solid_position_data=true)
Definition: temporary_stefan_boltzmann_elements.h:1247
oomph::StefanBoltzmannRadiationBase::sigma
double sigma()
Non-dim Stefan Boltzmann constant (switched off by default)
Definition: temporary_stefan_boltzmann_elements.h:1157
oomph::StefanBoltzmannRadiationBase::output_stefan_boltzmann_radiation_rays
void output_stefan_boltzmann_radiation_rays(std::ostream &outfile, const unsigned &integration_point=UINT_MAX)
Definition: temporary_stefan_boltzmann_elements.h:1487
oomph::StefanBoltzmannRadiationBase::sigma_pt
double *& sigma_pt()
Pointer to non-dim Stefan Boltzmann constant.
Definition: temporary_stefan_boltzmann_elements.h:1151
oomph::StefanBoltzmannRadiationBase::Theta_0_pt
double * Theta_0_pt
Pointer to non-dim zero centrigrade offset in Stefan Boltzmann law.
Definition: temporary_stefan_boltzmann_elements.h:1347
oomph::StefanBoltzmannRadiationBase::theta_0
double theta_0()
Definition: temporary_stefan_boltzmann_elements.h:1172
oomph::StefanBoltzmannRadiationBase::StefanBoltzmannRadiationBase
StefanBoltzmannRadiationBase()
Constructor.
Definition: temporary_stefan_boltzmann_elements.h:1137
oomph::StefanBoltzmannRadiationBase::contribution_to_stefan_boltzmann_radiation
double contribution_to_stefan_boltzmann_radiation(const Vector< double > &r_illuminated, const Vector< double > &n_illuminated, const Vector< unsigned > &visible_intpts_in_current_element)
Definition: temporary_stefan_boltzmann_elements.h:1199
oomph::StefanBoltzmannRadiationBase::wipe_stefan_boltzmann_illumination_info
void wipe_stefan_boltzmann_illumination_info()
Wipe illumination info.
Definition: temporary_stefan_boltzmann_elements.h:1225
oomph::StefanBoltzmannRadiationBase::theta_0_pt
double *& theta_0_pt()
Pointer to non-dim zero centrigrade offset in Stefan Boltzmann law.
Definition: temporary_stefan_boltzmann_elements.h:1154
oomph::StefanBoltzmannRadiationBase::output_stefan_boltzmann_radiation
void output_stefan_boltzmann_radiation(std::ostream &outfile)
Output Stefan Boltzmann radiation: x,y,in,out,n_x,n_y.
oomph::StefanBoltzmannRadiationBase::incoming_stefan_boltzmann_radiation
double incoming_stefan_boltzmann_radiation(const unsigned &ipt, const Vector< double > &r_illuminated, const Vector< double > &n_illuminated)
Definition: temporary_stefan_boltzmann_elements.h:1314
oomph::StefanBoltzmannRadiationBase::Sigma_pt
double * Sigma_pt
Pointer to non-dim Stefan Boltzmann constant.
Definition: temporary_stefan_boltzmann_elements.h:1344
oomph::StefanBoltzmannUnsteadyHeatFluxElement
Definition: temporary_stefan_boltzmann_elements.h:1572
oomph::StefanBoltzmannUnsteadyHeatFluxElement::StefanBoltzmannUnsteadyHeatFluxElement
StefanBoltzmannUnsteadyHeatFluxElement(FiniteElement *const &bulk_el_pt, const int &face_index)
Definition: temporary_stefan_boltzmann_elements.h:1770
oomph::StefanBoltzmannUnsteadyHeatFluxElement::StefanBoltzmannUnsteadyHeatFluxElement
StefanBoltzmannUnsteadyHeatFluxElement(const StefanBoltzmannUnsteadyHeatFluxElement &dummy)
Broken copy constructor.
Definition: temporary_stefan_boltzmann_elements.h:1582
oomph::StefanBoltzmannUnsteadyHeatFluxElement::output_stefan_boltzmann_radiation
void output_stefan_boltzmann_radiation(std::ostream &outfile)
Definition: temporary_stefan_boltzmann_elements.h:1613
oomph::StefanBoltzmannUnsteadyHeatFluxElement::get_flux
virtual void get_flux(const unsigned &ipt, const double &time, const Vector< double > &x, const Vector< double > &n, const double &u, double &flux)
Definition: temporary_stefan_boltzmann_elements.h:1721
oomph::StefanBoltzmannUnsteadyHeatFluxElement::set_integration_scheme
void set_integration_scheme(Integral *const &integral_pt)
Definition: temporary_stefan_boltzmann_elements.h:1590
oomph::StefanBoltzmannUnsteadyHeatFluxElement::zeta_nodal
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Definition: temporary_stefan_boltzmann_elements.h:1602
oomph::TemplateFreeUnsteadyHeatBaseFaceElement
Definition: heat_transfer_and_melt_elements.h:123
oomph::TemplateFreeUnsteadyHeatBaseFaceElement::Flux_fct_pt
UnsteadyHeatPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
Definition: heat_transfer_and_melt_elements.h:185
oomph::TemplateFreeUnsteadyHeatBaseFaceElement::U_index_ust_heat
unsigned U_index_ust_heat
Index at which temperature is stored.
Definition: heat_transfer_and_melt_elements.h:179
oomph::TemplateFreeUnsteadyHeatBaseFaceElement::Dim
unsigned Dim
The spatial dimension of the problem.
Definition: heat_transfer_and_melt_elements.h:182
oomph::TimeStepper::time_pt
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
oomph::Time::time
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
oomph::UnsteadyHeatBaseFaceElement
Definition: heat_transfer_and_melt_elements.h:1522
oomph::UnsteadyHeatBaseFaceElement::interpolated_u
void interpolated_u(const Vector< double > &s, double &u)
Temperature at local coordinate s.
Definition: heat_transfer_and_melt_elements.h:1669
oomph::Vector
Definition: oomph-lib/src/generic/Vector.h:58
min
#define min(a, b)
Definition: datatypes.h:22
max
#define max(a, b)
Definition: datatypes.h:23
s
RealScalar s
Definition: level1_cplx_impl.h:130
k
char char char int int * k
Definition: level2_impl.h:374
BiharmonicTestFunctions2::theta
double theta
Definition: two_d_biharmonic.cc:236
Eigen::bfloat16_impl::tanh
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 tanh(const bfloat16 &a)
Definition: BFloat16.h:639
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
GlobalFct::plot_it
void plot_it(const std::string filename)
Plot "landscape" of residuals (only for 2D problems!)
Definition: spring_contact.cc:211
MergeRestartFiles.filename
string filename
Definition: MergeRestartFiles.py:39
MergeRestartFiles.found
bool found
Definition: MergeRestartFiles.py:24
ProblemParameters::flux
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59
UniformPSDSelfTest.radius
radius
Definition: UniformPSDSelfTest.py:15
WallFunction::normal
void normal(const Vector< double > &x, Vector< double > &normal)
Definition: free_surface_rotation.cc:65
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
oomph::BrokenCopy::broken_copy
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Definition: oomph_utilities.cc:212
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
oomph::IntersectionChecker::intersects
bool intersects(const Vector< Vector< double > > &first_segment_vertex, const Vector< Vector< double > > &second_segment_vertex, const double &epsilon_parallel=1.0e-15)
Definition: temporary_stefan_boltzmann_elements.h:57
oomph::MathematicalConstants::Pi
const double Pi
50 digits from maple
Definition: oomph_utilities.h:157
oomph::SolarRadiationHelper::Zero_atmospheric_radiation_fct
void Zero_atmospheric_radiation_fct(const double &time, double &solar_flux_magnitude, Vector< double > &solar_flux_unit_vector, double &total_diffuse_radiation)
Default atmospheric radiation function in terms of time.
Definition: temporary_stefan_boltzmann_elements.h:123
oomph::StefanBoltzmannHelper::bin_helper
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)
Definition: temporary_stefan_boltzmann_elements.h:1861
oomph::StefanBoltzmannHelper::Nx_bin
unsigned Nx_bin
Number of bins in x direction.
Definition: temporary_stefan_boltzmann_elements.h:1833
oomph::StefanBoltzmannHelper::doc_sample_points
void doc_sample_points(std::ofstream &outfile, const Vector< FiniteElement * > &sb_face_element_pt)
Definition: temporary_stefan_boltzmann_elements.h:2124
oomph::StefanBoltzmannHelper::Dx
Vector< double > Dx(2, 0.0)
Increments in bin.
oomph::StefanBoltzmannHelper::Ny_bin
unsigned Ny_bin
Number of bins in y direction.
Definition: temporary_stefan_boltzmann_elements.h:1836
oomph::StefanBoltzmannHelper::Min_coord
Vector< double > Min_coord(2, DBL_MAX)
Min bin coords.
oomph::StefanBoltzmannHelper::Nsample
unsigned Nsample
Definition: temporary_stefan_boltzmann_elements.h:1840
oomph::StefanBoltzmannHelper::Max_coord
Vector< double > Max_coord(2,-DBL_MAX)
Max bin coords.
oomph::StefanBoltzmannHelper::doc_bins
void doc_bins(std::ofstream &bin_file)
Definition: temporary_stefan_boltzmann_elements.h:2093
oomph::StefanBoltzmannHelper::setup_stefan_boltzmann_visibility
void setup_stefan_boltzmann_visibility(const Vector< FiniteElement * > &sb_face_element_pt)
Definition: temporary_stefan_boltzmann_elements.h:2173
oomph::StefanBoltzmannHelper::Element_in_bin
Vector< Vector< std::set< FiniteElement * > > > Element_in_bin
Bin to store pointers to finite elements in.
Definition: temporary_stefan_boltzmann_elements.h:1821
oomph::TimingHelpers::timer
double timer()
returns the time in seconds after some point in past
Definition: oomph_utilities.cc:1295
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition: AnisotropicHookean.h:10
oomph::pause
void pause(std::string message)
Pause and display message.
Definition: oomph_utilities.cc:1265
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
plotDoE.x
list x
Definition: plotDoE.py:28
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2