#Example of geometric setup 3D frame structure
#unit system: force(N) - length(m) - mass(kg)

wipe; # clear memory of past model definitions
source utils_0.tcl

#.................................................................................................................................
#MODEL BUILDER

#3D model, 6dof
#model BasicBuilder -ndm $ndm <-ndf $ndf>
#(ndf and ndm are 1 based in OpenSees)
model BasicBuilder -ndm 3 -ndf 6

#.................................................................................................................................
#GEOMETRY

#structure
set L1	3.; 				#column height (parallel to Z axis)
set L2	[expr $L1*2]; 	#level 2
set L3	[expr $L1*3]; 	#level 3
set LX	4.; 				#beam length (parallel to X axis)
set LY	5.; 				#beam length (parallel to Y axis)

#create the nodes of the frame structure
#node $nodeTag (ndm $coords) <-mass (ndf $massValues)>
node 1		0.		0.		0.; #level0
node 2		$LX		0.		0.
node 3		$LX		$LY		0.
node 4		0.		$LY		0.
node 5		0.		0.		$L1; #level1
node 6		$LX		0.		$L1
node 7		$LX		$LY		$L1
node 8		0.		$LY		$L1
node 9		0.		0.		$L2; #level2
node 10		$LX		0.		$L2
node 11		$LX		$LY		$L2
node 12		0.		$LY		$L2
node 13		0.		0.		$L3; #level3
node 14		$LX		0.		$L3
node 15		$LX		$LY		$L3
node 16		0.		$LY		$L3

#create master nodes
node 17		2.		2.5		3.; #center level 0
node 18		2.		2.5		6.; #center level 1
node 19		2.		2.5		9.; #center level 2


#.................................................................................................................................
#MATERIALS
#tags
set ConcrTag 1; #unconfined cover concrete
set SteelTag 2; #reinforcement

#reinforced concrete
set fpc 	[expr -35.*$MPa]; 					# nominal concrete compressive strength (+tension, -compression)
set epsc0	-0.002;								# strain at maximum strength of unconfined concrete
set fpcu 	[expr 0.2*$fpc];					# ultimate stress
set epsU	-0.01;								# strain at ultimate stress
set lambda	0.1;								# ratio between unloading slope at $epsU and initial slope $Ec
set ft		[expr -0.14*$fpc];					# tensile strength +tension
set Ets 	[expr $ft/0.002];					# tension softening stiffness

#steel
set Fy 		[expr 300*$MPa];					# yield stress
set E 		[expr 210000.0*$MPa];				# young's modulus
set b 		0.05;								# strain-hardening ratio 
set R0		20;									# control the transition from elastic to plastic branches
set cR1		0.925
set cR2		0.15

#create the materials
uniaxialMaterial Concrete02	$ConcrTag $fpc $epsc0 $fpcu $epsU $lambda $ft $Ets;		# build cover concrete (unconfined)
uniaxialMaterial Steel02 $SteelTag $Fy $E $b $R0 $cR1 $cR2;						# build reinforcement material
puts "uniaxialMaterial Concrete02	$ConcrTag $fpc $epsc0 $fpcu $epsU $lambda $ft $Ets"
puts "uniaxialMaterial Steel02 $SteelTag $Fy $E $b $R0 $cR1 $cR2"


#.................................................................................................................................
#FIBER SECTIONS (symmetric)
#sections geometry and properties
set bCol	0.4
set hCol	0.4
set bBeam	0.3
set hBeam	0.5
set JCol 	[expr $bCol**4./12.]
set JBeam 	[expr $bBeam**4./12.]

set nu 0.30
set G [expr $E/(2.0*(1.0+$nu))]
set GJcol [expr $G*$JCol]
set GJbeam [expr $G*$JBeam]

set ColSecTag 		1;							# assign a tag number to the column section
set BeamSecTag 		2;							# assign a tag number to the column section

set diam			0.022;						# bars diameters
set cover 			[expr 0.03 + $diam/2.0];						# dimension of the concrete cover for both sections
set numBarsCol		5;							# total number of longitudinal-reinforcement bars in column. (symmetric top & bot)
set numBarsBeam		4;							# totalnumber of longitudinal-reinforcement bars in beam. (symmetric top & bot)
set barAreaCol		[expr 3.14*($diam/2)**2];	# area of longitudinal-reinforcement bars in columns
set barAreaBeam		[expr 3.14*($diam/2)**2];	# area of longitudinal-reinforcement bars in beams

set coverColY		[expr $bCol/2];				# distance from the column section z-axis to the edge of the cover concrete
set coverColZ		[expr $hCol/2];				# distance from the column section y-axis to the edge of the cover concrete
set coverBeamY		[expr $bBeam/2];			# distance from the beam section z-axis to the edge of the cover concrete
set coverBeamZ		[expr $hBeam/2];			# distance from the beam section y-axis to the edge of the cover concrete
set coreColY		[expr $coverColY-$cover];	# internal distance from column section z-axis to the top edge of confined concrete
set coreColZ		[expr $coverColZ-$cover];	# internal distance from column section y-axis to the top edge of confined concrete
set coreBeamY		[expr $coverBeamY-$cover];	# internal distance from beam section z-axis to the top edge of confined concrete
set coreBeamZ		[expr $coverBeamZ-$cover];	# internal distance from beam section y-axis to the top edge of confined concrete

set nfColY 			15;							# number of fibers for concrete column in y-direction
set nfCOlZ 			15;							# number of fibers for concrete column in z-direction
set nfBeamY 		9;							# number of fibers for concrete beam in y-direction
set nfBeamZ 		15;							# number of fibers for concrete beam in z-direction

#create the fiber section for columns
section Fiber $ColSecTag -GJ $GJcol {
	#concrete patch
	patch quadr $ConcrTag $nfCOlZ $nfColY -$coverColY $coverColZ -$coverColY -$coverColZ $coverColY -$coverColZ $coverColY $coverColZ
	# top layer reinfocement
	layer straight $SteelTag $numBarsCol $barAreaCol -$coreColY -$coreColZ $coreColY -$coreColZ
	# bottom layer reinforcement
	layer straight $SteelTag $numBarsCol $barAreaCol -$coreColY $coreColZ  $coreColY $coreColZ
}

#create the fiber section for beams
section Fiber $BeamSecTag -GJ $GJbeam {
	#concrete patch
	patch quadr $ConcrTag $nfBeamZ $nfBeamY -$coverBeamY $coverBeamZ -$coverBeamY -$coverBeamZ $coverBeamY -$coverBeamZ $coverBeamY $coverBeamZ
	# top layer reinfocement
	layer straight $SteelTag $numBarsBeam $barAreaBeam -$coreBeamY -$coreBeamZ $coreBeamY -$coreBeamZ
	# bottom layer reinforcement
	layer straight $SteelTag $numBarsBeam $barAreaBeam -$coreBeamY $coreBeamZ  $coreBeamY $coreBeamZ
}

#.................................................................................................................................
#ELEMENTS
#move from the basic CS to the global CS in 3D
#to match STKO default beam and column characteristic axis direction assignment use the geomTransf 1 for beams and 2 for columns
#geomTransf Linear $transfTag $vecxzX $vecxzY $vecxzZ <-jntOffset $dXi $dYi $dZi $dXj $dYj $dZj>
set ColTransfTag 	1
set BeamTransfTag 	2
geomTransf Linear $ColTransfTag 1 0 0
geomTransf Linear $BeamTransfTag 0 0 1

#integration points
set numIntgrPts 	5;							# number of integration points for force-based element

#2 loop structures to create elements of the frame
set k 				0;							#initialize a counter for the elements id

for {set i 1} {$i < 13} {incr i} {

		set node_bot $i; #bottom node of the column
		set node_top [expr $i+4]; #top node of the column
		
		#element forceBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag
		element forceBeamColumn $k $node_bot $node_top $numIntgrPts $ColSecTag $ColTransfTag
		incr k
}

set floor_id {0 1 2 3 0}
for {set i 1} {$i < 4} {incr i} {
		
		set index [expr $i*4 + 1]
		
		for {set j 0} {$j < 4} {incr j} {
			
			set n1 [expr $index+[lindex $floor_id $j]]
			set n2 [expr $index+[lindex $floor_id [expr $j+1]]]
			
			#element forceBeamColumn $eleTag $iNode $jNode $numIntgrPts $secTag $transfTag
			element forceBeamColumn $k $n1 $n2 $numIntgrPts $BeamSecTag $BeamTransfTag
			incr k
	}
}

#.................................................................................................................................
#DEFINITIONS
#create a timeseries
# timeSeries Linear
timeSeries Linear 1

#.................................................................................................................................
#CONDITIONS
#single-point constraints
#fix $nodeTag (ndf $constrValues)   
fix 1 1 1 1 1 1 1
fix 2 1 1 1 1 1 1
fix 3 1 1 1 1 1 1
fix 4 1 1 1 1 1 1
fix 17 0 0 1 1 1 0
fix 18 0 0 1 1 1 0
fix 19 0 0 1 1 1 0

#multi-point constraints
#rigidDiaphragm $perpDirn $rNodeTag $cNodeTag1 $cNodeTag2 ...
#$perpDirn direction perpendicular to the rigid plane (plane xy = 3, plane xz = 2, plane zy = 1)
#rNodeTag restrained node tag
#cNodeTags constrained node tags 
rigidDiaphragm 3 17 5 6 7 8
rigidDiaphragm 3 18 9 10 11 12
rigidDiaphragm 3 19 13 14 15 16

#loads
set p -2000.0
set n 0

pattern Plain 1 1 {
	for {set i 12} {$i < 24} {incr i} {
		eleLoad -ele $i -type -beamUniform 0.0 -2000.0 0.0
		puts "eleLoad -ele $i -type -beamUniform 0.0 -2000.0 0.0"
	}
}

puts "---------------------------"

#.................................................................................................................................
#ANALYSIS

#add the mpco recorder to compare the results with the STKO example (ex3.scd)
recorder mpco "output/ex8tcl.mpco" \
-N "displacement" "rotation" "reactionForce" \
-E "localForce" "section.force" "section.deformation" "section.fiber.stress" "section.fiber.strain"


#system systemType?
system UmfPack

#constraints constraintType?
constraints Penalty 1e+18 1e+18

#numberer numbererType?
numberer RCM

#test testType? (NormDispIncr $tol $iter)
test NormDispIncr 1e-06 10

#algorithm algorithmType? (Newton)
algorithm Newton

#integrator LoadControl $lambda
integrator LoadControl 0.1

#analysis analysisType?
analysis Static

#raise an error if the analysis does not work
#analyze $numIncr 
set ok [analyze 10]
if {$ok != 0} {
	error "Analysis failed!"
}

#print analysis done
puts "Static analysis done!"

#if you don't put this OpenSees will not save your data on the out put files until you exit it
wipe