#  two-beam 3D model: 
# problem: the displacement doesn't change under the heat
wipe;

model BasicBuilder -ndm 3 -ndf 6;

source Wsection.tcl; #PK Wsection includes fibersecThermal


source DisplayPlane.tcl
source DisplayModel2D.tcl
source DisplayModel3D.tcl
#Define nodes and coordinates
node 1 0. 0. 0.;
node 2 1000. 0. 0.;
node 3 2000. 0. 0.;
node 4 3000. 0. 0.;
node 5 4000. 0. 0.;
node 6 5000. 0. 0.;
node 7 6000. 0. 0.;

node 11 6000. 0. 1000.;
node 12 6000. 0. 2000.;
node 13 6000. 0. 3000.;
node 14 6000. 0. 4000.;
node 15 6000. 0. 5000.;
node 16 6000. 0. 6000.;


#all fixed 

fix 1 1 1 1 1 1 1;
fix 16 1 1 1 1 1 1; 



# =======  material and section
set secID 1
set matID 1
uniaxialMaterial Steel01Thermal 1 308 2.1e5 0.01; 

# beam sections: 
set d 355.;	# depth
set bf 171.5;	# flange width
set tf 11.5;	# flange thickness
set tw 7.4;	    # web thickness
set nfdw 16;		# number of fibers along dw
set nftw 2;		# number of fibers along tw
set nfbf 16;		# number of fibers along bf
set nftf 4;			# number of fibers along tf


#section
	set dw [expr $d - 2. * $tf]
	set y1 [expr -$d/2.]
	set y2 [expr -$dw/2.]
	set y3 [expr $dw/2.]
	set y4 [expr $d/2.]
	set z1 [expr -$bf/2.]
	set z2 [expr -$tw/2.]
	set z3 [expr $tw/2.]
	set z4 [expr $bf/2.]

	section FiberThermal $secID {
	patch quadr $matID $nfbf $nftf $y1 $z4 $y1 $z1 $y2 $z1 $y2 $z4
	patch quadr $matID $nftw $nfdw $y2 $z3 $y2 $z2 $y3 $z2 $y3 $z3
	patch quadr $matID $nfbf $nftf $y3 $z4 $y3 $z1 $y4 $z1 $y4 $z4
	}




      
# geomTransf Corotational 1;




#choose the section for the beam eles



set section 22
set Ubig 9.E15;
# set Ubig 9.E1;
set SecTagTorsion 3
uniaxialMaterial Elastic $SecTagTorsion $Ubig
section Aggregator $section $SecTagTorsion T -section $secID;


# ===================== element ==============================
set IDBeamXTransf 1;
set IDBeamZTransf 2;
set IDColTransf 3;
geomTransf Corotational $IDBeamXTransf 0 0 1 ;
geomTransf Corotational $IDBeamZTransf -1 0 0 ;  #1 0 0
geomTransf Corotational $IDColTransf 0 0 1 ;

element dispBeamColumnThermal 1 1 2 5 $section 1;
element dispBeamColumnThermal 2 2 3 5 $section 1;
element dispBeamColumnThermal 3 3 4 5 $section 1;
element dispBeamColumnThermal 4 4 5 5 $section 1;
element dispBeamColumnThermal 5 5 6 5 $section 1;
element dispBeamColumnThermal 6 6 7 5 $section 1;

 element dispBeamColumnThermal 11 7 11 5 $section 2;
element dispBeamColumnThermal 12 11 12 5 $section 2;
element dispBeamColumnThermal 13 12 13 5 $section 2;
element dispBeamColumnThermal 14 13 14 5 $section 2;
element dispBeamColumnThermal 15 14 15 5 $section 2;
element dispBeamColumnThermal 16 15 16 5 $section 2;



recorder Node -file 7.out -time -node 7 -dof 1 2 3 4 5 6 disp;


# Define DISPLAY -------------------------------------------------------------
	set  xPixels 700;	# height of graphical window in pixels
	set  yPixels 400;	# height of graphical window in pixels
	set  xLoc1 10;	# horizontal location of graphical window (0=upper left-most corner)
	set  yLoc1 400;	# vertical location of graphical window (0=upper left-most corner)
	set ViewScale 0.0000001;	# scaling factor for viewing deformed shape, it depends on the dimensions of the model
	DisplayModel3D DeformedShape $ViewScale $xLoc1 $yLoc1 $xPixels $yPixels 0

pattern Plain 1 Linear {

#CREATE UNIFORM LOADS FOR BEAMS

load 7 0.  -100. 0. 0. 0. 0. 
}

	constraints Plain;					# how it handles boundary conditions
	##numberer Plain;	
	numberer RCM;					# renumber dof's to minimize band-width (optimization)
	#system BandGeneral;					# how to store and solve the system of equations in the analysis (large model: try UmfPack)
	system UmfPack
	test NormUnbalance 1.0e-5 12;		# tolerance, max iterations
	algorithm Newton;					# use Newton's solution algorithm: updates tangent stiffness at every iteration
	
	integrator LoadControl 0.1;				# determine the next time step for an analysis, # apply gravity in 10 steps
	analysis Static;					# define type of analysis static or transient
	analyze 10;
	
	print domain.out
	loadConst -time 0.0	;

	
	
	
	
	
	
	
	#######FIRE PART########
	puts "Fire";
	
	
	set HalfD [expr $d/2.0];
	set HalfB [expr $bf/2.0];
	#create fire pattern
pattern Plain 2 Linear {
  # not sure how 3dthermalAction is exactly set: I try all the following, but fail
	# eleLoad -range 1 6 -type -beam3dThermal -source Tdata.dat -$HalfD  $HalfD -$HalfB $HalfB; 
	# eleLoad -range 11 16 -type -beam3dThermal -source Tdata.dat -$HalfD  $HalfD -$HalfB $HalfB; 
	
	# eleLoad -range 1 6 -type -beam3dThermal 1000 -$HalfD  1000 $HalfD 
	# eleLoad -range 11 16 -type -beam3dThermal  1000 -$HalfD 1000 $HalfD 
	
	eleLoad -range 1 6 -type -beam3dThermalAction 1000 -$HalfD  1000 $HalfD 
	eleLoad -range 11 16 -type -beam3dThermalAction  1000 -$HalfD 1000 $HalfD 
	
}

	#constraints Transformation
	#constraints Penalty 1e10 1e10
	constraints Plain;					# how it handles boundary conditions
	numberer Plain;						# renumber dof's to minimize band-width (optimization)
	#system BandGeneral;					# how to store and solve the system of equations in the analysis (large model: try UmfPack)
	system UmfPack
	test NormUnbalance 1.0e-4 80;		# tolerance, max iterations
	#test NormDispIncr 1e-6 10
	#test EnergyIncr 1.0e-8 50
	algorithm Newton;					# use Newton's solution algorithm: updates tangent stiffness at every iteration
	integrator LoadControl 0.01;
    analysis Static;
    analyze 100; #100steps to apply thermal action
	



	

set Dstep 0.0;
set ok 0
set Tol 1.0e-8;
set maxNumIter 100;
set printFlag 0;
set NewmarkGamma 0.5;
set NewmarkBeta 0.25;
set TestType NormDispIncr; # options: NormDispIncr EnergyIncr
set IDctrlNode 7;
set IDctrlDOF 2;

set Dmax 200;

variable TolStatic 1.e-8;                        # Convergence Test: tolerance
variable maxNumIterStatic 6;                # Convergence Test: maximum number of iterations that will be performed before "failure to converge" is returned
variable printFlagStatic 0;                # Convergence Test: flag used to print information on convergence (optional)        # 1: print information on each step; 
variable testTypeStatic EnergyIncr ;	# Convergence-test type
test $testTypeStatic $TolStatic $maxNumIterStatic $printFlagStatic;

variable algorithmTypeStatic Newton
algorithm $algorithmTypeStatic;     

while {$Dstep <= 1.0 && $ok == 0} {	  # convengent and go on
	set controlDisp [nodeDisp $IDctrlNode $IDctrlDOF ];
	set currentTime [getTime];
	puts "controlDisp:    ******************************************************************************************************************************            $currentTime     $controlDisp  $Dmax";
	set Dstep [expr $controlDisp/$Dmax];
	set ok [analyze 1];
	# puts "ok = $ok";
	set StaticConvProc naive
	if { $ok !=0} {

			if {$ok != 0} {  
					# if analysis fails, we try some other stuff
					# performance is slower inside this loop	global maxNumIterStatic;	    # max no. of iterations performed before "failure to converge" is ret'd
					puts "Trying Newton with Initial Tangent .."
					test NormDispIncr   $Tol 2000 0
					algorithm Newton -initial
					set ok [analyze 1]
					test $testTypeStatic $TolStatic  $maxNumIterStatic   0
					algorithm $algorithmTypeStatic
			}
			if {$ok != 0} {
							puts "Trying Broyden .."
							algorithm Broyden 8
							set ok [analyze 1 ]
							algorithm $algorithmTypeStatic
						}
			if {$ok != 0} {
							puts "Trying NewtonWithLineSearch .."
							algorithm NewtonLineSearch 0.8 
							set ok [analyze 1]
							algorithm $algorithmTypeStatic
			}

	}

}