{  GOS_FlexPDE_2D_Beam_01.pde
}

title 'Beam 3D heat conduction'

select
    regrid=off { use fixed grid }
    ngrid=30     { smaller grid for quicker run }

coordinates
    cartesian2

variables
    Tpc (threshold=0.01)   { the temperature variable, with approximate size  }

global variables
	tpcmax
	tpcmax_t
	tpcmax_x
	tpcmax_y

time 0 to 100000 by 60 { establish time range and initial timestep }

definitions

    timeMax = 50400
    Ta = 20 + 273
    bottomTemp = Ta
    surfaceTemp = Ta
    length = 1					{ length }
    wide = 1.2					{ width }
    thick = 0.2					{ thickness }
    beamTopY = thick/2
    beamTopLowerY = beamTopY - 0.05

    K = 1.3  						{ Thermal conductivity  }
    H = 6
	Cp = 880 
    m = 2400 ! this is really rho, the density of the material !
    C = Cp * m

 Q = 20 * Cp * Sin(2*pi* t / 36000)  * upulse(t , t - 18000)   { Thermal source }

! Cannot get a max along this X axis
!TpcSurfaceFunc = Passive( IF ((y <= beamTopY) AND (y >= beamTopLowerY)) THEN Tpc ELSE TpcSurfaceMax )
!TpcSurfaceFunc = IF ((y <= beamTopY) AND (y >= beamTopLowerY)) THEN Tpc ELSE TpcSurfaceMax 
! Can get a max of Tpc, but this is not what I want
! TpcSurfaceFunc(y1) = Tpc
!TpcSurfaceMaxFunc = Passive( max(TpcSurface, TpcSurfaceMax) )
!TpcSurfaceMaxFunc = max(TpcSurface, TpcSurfaceMax)

tpcmax1 = time_max(globalmax(tpc))
tpgmax = globalmax(tpc)
tpcmaxfun = if(tpgmax > tpcmax) then tpgmax else tpcmax;
tpcmax_tfun = if(tpgmax > tpcmax) then t else tpcmax_t;
tpcmax_xfun = if(tpgmax > tpcmax) then globalmax_x(tpc) else tpcmax_x;
tpcmax_yfun = if(tpgmax > tpcmax) then globalmax_y(tpc) else tpcmax_y;

Tp = Tpc + 273


initial values	
    Tpc = 20

equations
    Tpc: div(k*grad(Tpc)) + Q = C * dt(Tpc) 	{ the heat equation }  
	tpcmax: tpcmax = tpcmaxfun;
	tpcmax_t: tpcmax_t = tpcmax_tfun
	tpcmax_x: tpcmax_x = tpcmax_xfun;
	tpcmax_y: tpcmax_y = tpcmax_yfun;

 boundaries

    Region 'BeamRegion'    							{ Define full domain boundary }
       start "Beam" (-wide/2, -thick/2)
		 	natural(Tpc) = -H * 0.5 *(Tp - Ta)  		{ Set convection value for Bottom }
         line to  (wide/2, -thick/2)  			
		 	natural(Tpc) = -H * 1 * (Tp - Ta)  		{ Set convection value for Right }
         line to  (wide/2, thick/2)
		 	natural(Tpc) = -H * 1 *(Tp - Ta)  		{ Set convection value for Top }
         line to (-wide/2 , thick/2)
		 	natural(Tpc) = -H * 1 * (Tp - Ta)  		{ Set convection value for Left }
         line to (-wide/2, -thick/2)


monitors
  for cycle=1
    contour(Tpc) on "BeamRegion" as "Beam Temp"  painted  range=(0,tpcmax)

plots
  for t = timeMax 
    contour(Tpc) on "BeamRegion" as "Beam Temp at Tmax" painted  range=(0,tpcmax)

    elevation(Tpc) from (-wide/2, thick/2) to (wide/2, thick/2) as "Top Surface Temp"


histories
    history(Tpc) at 
                    (0, thick/2)  ! range=(0,tmax)
                    (0, 0)  ! range=(0,tmax)
                    (0, -thick/2)   range=(0,tpcmax)
		report(tpcmax) report(tpcmax_t) report(tpcmax_x) report(tpcmax_y)

    history(tpgmax,tpcmax,tpcmax1) 
		report(tpcmax) report(tpcmax_t) report(tpcmax_x) report(tpcmax_y)


end