Three-Dimensional rigid Frame subjected to bi-directional earthquake ground motion

Figure 1: Three-Dimensional Rigid Frame & horizontal displacements at the diaphragm master nodes

Time History Analysis


//# Units: kip, in, sec
//# ----------------------------
//# Start of model generation
//# ----------------------------
//# Create ModelBuilder with 3 dimensions and 6 DOF/node
//model BasicBuilder -ndm 3 - ndf 6
var theModel = new OpenSees.AnalysisModelWrapper();


var theDomain = new OpenSees.Components.DomainWrapper();

//# Define geometry
//# ---------------

//# Set parameters for model geometry
//            set h 144.0;       # Story height
//set by 240.0;      # Bay width in Y-direction
//set bx 240.0;      # Bay width in X-direction


var h = 144.0;       //Story height
var by = 240.0;      //Bay width in Y-direction
var bx = 240.0;      //Bay width in X-direction


//# Create nodes
//# tag             X             Y           Z 
//        node  1[expr -$bx / 2][expr  $by / 2]           0
//node  2[expr  $bx / 2][expr  $by / 2]           0
//node  3[expr  $bx / 2][expr -$by / 2]           0
//node  4[expr -$bx / 2][expr -$by / 2]           0

//node  5[expr -$bx / 2][expr  $by / 2]          $h
//node  6[expr  $bx / 2][expr  $by / 2]          $h
//node  7[expr  $bx / 2][expr -$by / 2]          $h
//node  8[expr -$bx / 2][expr -$by / 2]          $h

//node 10[expr -$bx / 2][expr  $by / 2][expr 2 *$h]
//node 11[expr  $bx / 2][expr  $by / 2][expr 2 *$h]
//node 12[expr  $bx / 2][expr -$by / 2][expr 2 *$h]
//node 13[expr -$bx / 2][expr -$by / 2][expr 2 *$h]

//node 15[expr -$bx / 2][expr  $by / 2][expr 3 *$h]
//node 16[expr  $bx / 2][expr  $by / 2][expr 3 *$h]
//node 17[expr  $bx / 2][expr -$by / 2][expr 3 *$h]
//node 18[expr -$bx / 2][expr -$by / 2][expr 3 *$h]

theDomain.AddNode(new OpenSees.Components.NodeWrapper(1, 6, -bx / 2, by / 2, 0));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(2, 6, bx / 2, by / 2, 0));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(3, 6, bx / 2, -by / 2, 0));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(4, 6, -bx / 2, -by / 2, 0));

theDomain.AddNode(new OpenSees.Components.NodeWrapper(5, 6, -bx / 2, by / 2, h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(6, 6, bx / 2, by / 2, h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(7, 6, bx / 2, -by / 2, h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(8, 6, -bx / 2, -by / 2, h));

theDomain.AddNode(new OpenSees.Components.NodeWrapper(10, 6, -bx / 2, by / 2, 2 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(11, 6, bx / 2, by / 2, 2 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(12, 6, bx / 2, -by / 2, 2 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(13, 6, -bx / 2, -by / 2, 2 * h));

theDomain.AddNode(new OpenSees.Components.NodeWrapper(15, 6, -bx / 2, by / 2, 3 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(16, 6, bx / 2, by / 2, 3 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(17, 6, bx / 2, -by / 2, 3 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(18, 6, -bx / 2, -by / 2, 3 * h));

//# Master nodes for rigid diaphragm
//# tag X Y          Z 
//node  9  0 0         $h
//node 14  0 0[expr 2 *$h]
//node 19  0 0[expr 3 *$h]
theDomain.AddNode(new OpenSees.Components.NodeWrapper(9, 6, 0, 0, 1 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(14, 6, 0, 0, 2 * h));
theDomain.AddNode(new OpenSees.Components.NodeWrapper(19, 6, 0, 0, 3 * h));

//# Set base constraints
//# tag DX DY DZ RX RY RZ
//            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

foreach (var nodeId in new int[] { 1, 2, 3, 4 })
{
	for (var i = 0; i < 6; i++)
	{
		theDomain.AddSP_Constraint(new OpenSees.Components.Constraints.SP_ConstraintWrapper(nodeId, i, 0.0, true));
	}
}

//# Define rigid diaphragm multi-point constraints
//# normalDir  master     slaves
//rigidDiaphragm     3          9     5  6  7  8
//rigidDiaphragm     3         14    10 11 12 13
//rigidDiaphragm     3         19    15 16 17 18

theDomain.CreateRigidDiaphragm(9, new IDWrapper(new int[] { 5, 6, 7, 8 }), 2);
theDomain.CreateRigidDiaphragm(14, new IDWrapper(new int[] { 10, 11, 12, 13 }), 2);
theDomain.CreateRigidDiaphragm(19, new IDWrapper(new int[] { 15, 16, 17, 18 }), 2);

//# Constraints for rigid diaphragm master nodes
//# tag DX DY DZ RX RY RZ
//fix  9   0  0  1  1  1  0
//fix 14   0  0  1  1  1  0
//fix 19   0  0  1  1  1  0
foreach (var nodeId in new int[] { 9, 14, 19 })
{
	for (var i = 2; i < 5; i++)
	{
		theDomain.AddSP_Constraint(new OpenSees.Components.Constraints.SP_ConstraintWrapper(nodeId, i, 0.0, true));
	}
}

//# Define materials for nonlinear columns
//# --------------------------------------
//# CONCRETE
//# Core concrete (confined)
//# tag  f'c  epsc0  f'cu  epscu
//uniaxialMaterial Concrete01  1 - 5.0 - 0.005 - 3.5 - 0.02
var mat1 = new OpenSees.Materials.Uniaxials.Concrete01Wrapper(1, -5.0, -0.005, -3.5, -0.02);

//# Cover concrete (unconfined)
var fc = 4.0;
var mat2 = new OpenSees.Materials.Uniaxials.Concrete01Wrapper(2, -fc, -0.002, 0, -0.006);

//# STEEL
//# Reinforcing steel
//# tag fy   E     b
//uniaxialMaterial Steel01  3  60 30000 0.02
var mat3 = new OpenSees.Materials.Uniaxials.Steel01Wrapper(3, 60, 30000, 0.02);

var colh = 18.0;
var colb = 18.0;

//# Concrete elastic stiffness
//set E[expr 57000.0 * sqrt($fc * 1000) / 1000];
var Ec = 57000.0 * Math.Sqrt(fc * 1000) / 1000;
var GJ = 1.0e10;

var GJMat = new OpenSees.Materials.Uniaxials.ElasticMaterialWrapper(10, GJ, 0);

OpenSees.Materials.Sections.SectionForceDeformationWrapper colSection = null;
if (!elastic)
{
	colSection = CreateRCSection(1, colh, colb, 2.5,
	new KeyValuePair<int, Materials.Uniaxials.UniaxialMaterialWrapper>(1, mat1),
	new KeyValuePair<int, Materials.Uniaxials.UniaxialMaterialWrapper>(2, mat2),
	new KeyValuePair<int, Materials.Uniaxials.UniaxialMaterialWrapper>(3, mat3),
	3, 0.79, 8, 8, 10, 10, GJMat);
}
else

{
	colSection = new OpenSees.Materials.Sections.ElasticSection3dWrapper(4, Ec, colb * colh, 1 / 12.0 * colb * Math.Pow(colh, 3), 1 / 12.0 * colh * Math.Pow(colb, 3), GJ, 1);
}


//#define geometric transformation: performs a linear geometric transformation of beam stiffness and resisting force from the basic system to the global-coordinate system
//set ColTransfTag 1; 			# associate a tag to column transformation
//set BeamTransfTag 2; 			# associate a tag to beam transformation (good practice to keep col and beam separate)
//set ColTransfType Linear;			# options, Linear PDelta Corotational 
//geomTransf $ColTransfType $ColTransfTag; 	# only columns can have PDelta effects (gravity effects)
//geomTransf Linear $BeamTransfTag;
var pdelta = false;
Elements.CrdTransfs.CrdTransfWrapper colTransf = null;
if (!pdelta)
	colTransf = new Elements.CrdTransfs.LinearCrdTransf3dWrapper(new VectorWrapper(new double[] { 1, 0, 0 })); 
else
	colTransf = new Elements.CrdTransfs.PDeltaCrdTransf3dWrapper(new VectorWrapper(new double[] { 1, 0, 0 }));


//# Number of column integration points (sections)
//set np 4	
//# Create the nonlinear column elements
//# tag ndI ndJ nPts   secID  transf
//element nonlinearBeamColumn  1   1   5   $np  $colSec    1
//element nonlinearBeamColumn  2   2   6   $np  $colSec    1
//element nonlinearBeamColumn  3   3   7   $np  $colSec    1
//element nonlinearBeamColumn  4   4   8   $np  $colSec    1

//element nonlinearBeamColumn  5   5  10   $np  $colSec    1
//element nonlinearBeamColumn  6   6  11   $np  $colSec    1
//element nonlinearBeamColumn  7   7  12   $np  $colSec    1
//element nonlinearBeamColumn  8   8  13   $np  $colSec    1

//element nonlinearBeamColumn  9  10  15   $np  $colSec    1
//element nonlinearBeamColumn 10  11  16   $np  $colSec    1
//element nonlinearBeamColumn 11  12  17   $np  $colSec    1
//element nonlinearBeamColumn 12  13  18   $np  $colSec    1

var np = 4;

var integrationRule = new Elements.BeamIntegrations.BeamIntegrationWrapper(Elements.BeamIntegrations.BeamIntegrationType.Lobatto);
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(1, 1, 5, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(2, 2, 6, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(3, 3, 7, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(4, 4, 8, np, colSection, integrationRule, colTransf));

theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(5, 5, 10, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(6, 6, 11, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(7, 7, 12, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(8, 8, 13, np, colSection, integrationRule, colTransf));

theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(9, 10, 15, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(10, 11, 16, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(11, 12, 17, np, colSection, integrationRule, colTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(12, 13, 18, np, colSection, integrationRule, colTransf));

//# Define beam elements
//# --------------------
var Abeam = 18 * 24;
var Ibeamzz = 0.5 * 1.0 / 12 * 18.0 * Math.Pow(24, 3);
var Ibeamyy = 0.5 * 1.0 / 12 * 24.0 * Math.Pow(18, 3);

//# Define elastic section for beams
//# tag  E    A      Iz       Iy     G   J
//section Elastic  3  $E $Abeam $Ibeamzz $Ibeamyy $GJ 1.0
var beamSection = new OpenSees.Materials.Sections.ElasticSection3dWrapper(3, Ec, Abeam, Ibeamzz, Ibeamyy, GJ, 1);

//# Geometric transformation for beams
//# tag  vecxz
//geomTransf Linear 2   1 1 0
Elements.CrdTransfs.LinearCrdTransf3dWrapper beamTransf = new Elements.CrdTransfs.LinearCrdTransf3dWrapper(new VectorWrapper(new double[] { 1, 1, 0 }));

var beam_np = 5;

//# Create the beam elements
//# tag ndI ndJ nPts    secID   transf
//element nonlinearBeamColumn  13  5   6   $np  $beamSec     2
//element nonlinearBeamColumn  14  6   7   $np  $beamSec     2
//element nonlinearBeamColumn  15  7   8   $np  $beamSec     2
//element nonlinearBeamColumn  16  8   5   $np  $beamSec     2

//element nonlinearBeamColumn  17 10  11   $np  $beamSec     2
//element nonlinearBeamColumn  18 11  12   $np  $beamSec     2
//element nonlinearBeamColumn  19 12  13   $np  $beamSec     2
//element nonlinearBeamColumn  20 13  10   $np  $beamSec     2

//element nonlinearBeamColumn  21 15  16   $np  $beamSec     2
//element nonlinearBeamColumn  22 16  17   $np  $beamSec     2
//element nonlinearBeamColumn  23 17  18   $np  $beamSec     2
//element nonlinearBeamColumn  24 18  15   $np  $beamSec     2

theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(13, 5, 6, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(14, 6, 7, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(15, 7, 8, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(16, 8, 5, beam_np, beamSection, integrationRule, beamTransf));

theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(17, 10, 11, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(18, 11, 12, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(19, 12, 13, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(20, 13, 10, beam_np, beamSection, integrationRule, beamTransf));

theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(21, 15, 16, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(22, 16, 17, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(23, 17, 18, beam_np, beamSection, integrationRule, beamTransf));
theDomain.AddElement(new Elements.ForceBeamColumn3dWrapper(24, 18, 15, beam_np, beamSection, integrationRule, beamTransf));

//# Define gravity loads
//# --------------------
//# Gravity load applied at each corner node
//# 10% of column capacity
var p = 0.1 * fc * colh * colh;
var g = 386.4;
var mass = (4 * p) / g;

//# Rotary inertia of floor about master node
var inertia = mass * (bx * bx + by * by) / 12.0;

//# Set mass at the master nodes
//# tag MX MY MZ RX RY RZ
//mass  9  $m $m  0  0  0 $i
//mass 14  $m $m  0  0  0 $i
//mass 19  $m $m  0  0  0 $i

var massMatrix = new OpenSees.MatrixWrapper(6, 6);
massMatrix.Zero();
massMatrix[0, 0] = massMatrix[1, 1] = mass;
massMatrix[5, 5] = inertia;
theDomain.SetMass(new int[] { 9, 14, 19 }, massMatrix);


//# Define gravity loads
//pattern Plain 1 Constant {
//    foreach node { 5 6 7 8  10 11 12 13  15 16 17 18}
//    {
//        load $node 0.0 0.0 -$p 0.0 0.0 0.0
//    }
//}
var theSeries = new OpenSees.Components.Timeseries.ConstantSeriesWrapper(1.0);
var theLoadPattern = new OpenSees.Components.LoadPatterns.LoadPatternWrapper(1);
theLoadPattern.SetTimeSeries(theSeries);
theDomain.AddLoadPattern(theLoadPattern);
var nltag = 1;
foreach (var node in new int[] { 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 18 })
{
	theDomain.AddNodalLoad(new Components.Loads.NodalLoadWrapper(nltag++, node, new VectorWrapper(new double[] { 0, 0, -p, 0, 0, 0 }), false), 1);
}


theDomain.SetRayleighDampingFactors(0, 0, 0, 0.0018);

//# Define earthquake excitation
//# ----------------------------
//# Set up the acceleration records for Tabas fault normal and fault parallel
//set tabasFN "Path -filePath tabasFN.txt -dt 0.02 -factor $g"
//set tabasFP "Path -filePath tabasFP.txt -dt 0.02 -factor $g"
//# Define the excitation using the Tabas ground motion records
//# tag dir         accel series args
//pattern UniformExcitation  2   1 - accel      $tabasFN
//pattern UniformExcitation  3   2 - accel      $tabasFP

{
	var dt = 0.02;
	var factor = g;
	var accelVector = GetFileDataSeries([email protected]"{Environment.CurrentDirectory}\tabasFN.txt");
	var accelSeries = new OpenSees.Components.Timeseries.PathSeriesWrapper(accelVector, dt, factor, false, false, 0);
	var groundMotion = new OpenSees.Components.GroundMotions.GroundMotionWrapper(null, null, accelSeries, null, dt, 1);
	var uniformExcitationPattern = new OpenSees.Components.LoadPatterns.UniformExcitationWrapper(2, groundMotion, 0, 0, 1);
	theDomain.AddLoadPattern(uniformExcitationPattern);
}

{
	var dt = 0.02;
	var factor = g;
	var accelVector = GetFileDataSeries([email protected]"{Environment.CurrentDirectory}\tabasFP.txt");
	var accelSeries = new OpenSees.Components.Timeseries.PathSeriesWrapper(accelVector, dt, factor, false, false, 0);
	var groundMotion = new OpenSees.Components.GroundMotions.GroundMotionWrapper(null, null, accelSeries, null, dt, 1);
	var uniformExcitationPattern = new OpenSees.Components.LoadPatterns.UniformExcitationWrapper(3, groundMotion, 1, 0, 1);
	theDomain.AddLoadPattern(uniformExcitationPattern);
}

//# Create the convergence test
//# tol   maxIter  printFlag
//test EnergyIncr 1.0e-8   20         3
var tol = 1.0e-8;
var theTest = new OpenSees.ConvergenceTests.CTestEnergyIncrWrapper(tol, 20, 3, 2);
var theSolnAlgo = new OpenSees.Algorithms.NewtonRaphsonWrapper();
var theHandler = new OpenSees.Handlers.TransformationConstraintHandlerWrapper();

var theSolver = new OpenSees.Systems.Linears.SuperLUWrapper();
var theSOE = new OpenSees.Systems.Linears.SparseGenColLinSOEWrapper(theSolver);
var theNumberer = new OpenSees.Numberers.PlainNumbererWrapper();

var theIntegrator = new OpenSees.Integrators.Transient.NewmarkWrapper(0.5, 0.25);
var theTransientAnalysis = new OpenSees.Analysis.DirectIntegrationAnalysisWrapper(
	  theDomain,
	  theHandler,
	  theNumberer,
	  theModel,
	  theSolnAlgo,
	  theSOE,
	  theIntegrator,
	  theTest
  );



// eigen analysis
var theEigenSOE = new OpenSees.Systems.Eigens.FullGenEigenSOEWrapper(new OpenSees.Systems.Eigens.FullGenEigenSolverWrapper(), theModel);
theTransientAnalysis.SetEigenSOE(theEigenSOE);
// puts "eigen values at start of transient: [eigen 2]"
var ret = theTransientAnalysis.Eigen(3, true, true);
if (ret == 0)
{
	var eigens = theDomain.GetEigenValues();
	Console.WriteLine($"eigen values at start of transient: {string.Join(" ", eigens)}");
}
else
{
	Console.WriteLine($"eigen values analysis failed, return");
	return;
}

var recdofs = new OpenSees.IDWrapper(new int[] { 0, 1 });

{
	var opsstream = new OpenSees.Handlers.DataFileStreamWrapper(Environment.CurrentDirectory + $"\\Node51.out");
	var recorder = new OpenSees.Recorders.NodeRecorderWrapper(
		recdofs,
		new OpenSees.IDWrapper(new int[] { 9, 14, 19 }),
		0,
		"disp",
		theDomain,
		opsstream
	);
	theDomain.AddRecorder(recorder);
}

//#define type of analysis static or transient
//analyze $NstepGravity;		# apply gravity
ret = theTransientAnalysis.Analyze(2000, 0.01);
if (ret != 0)
{
	Console.WriteLine("analysis FAILED");
	Environment.Exit(-1);
}
else
{
	Console.WriteLine("analysis SUCCESSFULL");
}

theTransientAnalysis.ClearAll();
return;

static OpenSees.VectorWrapper GetFileDataSeries(string filename, int skipline = 0, double factor = 1)
{
	var sr = new System.IO.StreamReader(filename);
	for (var i = 0; i < skipline; i++)
		sr.ReadLine();
	var data = new List();
	while (!sr.EndOfStream)
	{
		var lineData = sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries).Select(i => double.Parse(i) * factor);
		if (lineData.Count() > 0)
			data.AddRange(lineData);
	}
	return new OpenSees.VectorWrapper(data.ToArray());
}

public static OpenSees.Materials.Sections.SectionForceDeformationWrapper CreateRCSection(
	int id,
	double h,
	double b,
	double cover,
	KeyValuePair<int, Materials.Uniaxials.UniaxialMaterialWrapper> coreMat,
	KeyValuePair<int, Materials.Uniaxials.UniaxialMaterialWrapper> coverMat,
	KeyValuePair<int, Materials.Uniaxials.UniaxialMaterialWrapper> rebarMat,
	int numBars,
	double barArea,
	int nfCoreY, int nfCoreZ, int nfCoverY, int nfCoverZ, Materials.Uniaxials.UniaxialMaterialWrapper GJ
	)
{
	var coverY = h / 2.0;
	var coverZ = b / 2.0;
	var ncoverY = -coverY;
	var ncoverZ = -coverZ;

	var coreY = coverY - cover;
	var coreZ = coverZ - cover;
	var ncoreY = -coreY;
	var ncoreZ = -coreZ;

	//Define the fiber section
	//var spacingY = (coreY - ncoreY) / (numBars - 1);
	//var numYBars = numBars - 2;
	var repres = new OpenSees.Materials.Sections.Repres.FiberSectionReprWrapper(1,
	   new OpenSees.Materials.Sections.Repres.PatchWrapper[]
	   {
		   
			// Define the core patch
			new OpenSees.Materials.Sections.Repres.QuadPatchWrapper(coreMat.Key,nfCoreZ,nfCoreY, new OpenSees.MatrixWrapper(new double[,]{{ ncoreY, coreZ },{ ncoreY, ncoreZ },{ coreY, ncoreZ },{ coreY, coreZ } })),

			////// Define the four cover patches
			////// patch quadr $coverID 1 $nfCoverY $ncoverY $coverZ $ncoreY $coreZ $coreY $coreZ $coverY $coverZ
			new OpenSees.Materials.Sections.Repres.QuadPatchWrapper(coverMat.Key,1,nfCoverY,new OpenSees.MatrixWrapper(new double[,]{{ ncoverY, coverZ },{ ncoreY, coreZ },{ coreY, coreZ },{ coverY, coverZ } })),
			//////// patch quadr $coverID 1 $nfCoverY $ncoreY $ncoreZ $ncoverY $ncoverZ $coverY $ncoverZ $coreY $ncoreZ                                                                                                                                                                            
			new OpenSees.Materials.Sections.Repres.QuadPatchWrapper(coverMat.Key,1,nfCoverY,new OpenSees.MatrixWrapper(new double[,]{{ ncoreY, ncoreZ }, { ncoverY, ncoverZ }, {coverY, ncoverZ }, { coreY, ncoreZ } })),
			//////// patch quadr $coverID $nfCoverZ 1 $ncoverY $coverZ $ncoverY $ncoverZ $ncoreY $ncoreZ $ncoreY $coreZ                                                                                                                                                                                                                                                                          
			new OpenSees.Materials.Sections.Repres.QuadPatchWrapper(coverMat.Key,nfCoverZ,1,new OpenSees.MatrixWrapper(new double[,]{ { ncoverY, coverZ }, { ncoverY, ncoverZ }, {ncoreY, ncoreZ }, { ncoreY, coreZ } })),
			////////patch quadr $coverID $nfCoverZ 1 $coreY $coreZ $coreY $ncoreZ $coverY $ncoverZ $coverY $coverZ
			new OpenSees.Materials.Sections.Repres.QuadPatchWrapper(coverMat.Key,nfCoverZ,1,new OpenSees.MatrixWrapper(new double[,]{ {coreY, coreZ }, { coreY, ncoreZ }, {coverY, ncoverZ }, { coverY, coverZ } })),
	   },
	   new OpenSees.Materials.Sections.Repres.ReinfLayerWrapper[]
	   {
				//# Define the steel along constant values of y (in the z direction)
				// layer straight $steelID $numBars $barArea $ncoreY $coreZ $ncoreY $ncoreZ
				//new OpenSees.Materials.Sections.Repres.StraightReinfLayerWrapper(rebarMat.Key,numBars,barArea,new OpenSees.VectorWrapper(new double[]{ ncoreY, coreZ }), new OpenSees.VectorWrapper(new double[]{ ncoreY, ncoreZ})),
				//// layer straight $steelID $numBars $barArea $coreY $coreZ $coreY $ncoreZ
				//new OpenSees.Materials.Sections.Repres.StraightReinfLayerWrapper(rebarMat.Key,numBars,barArea,new OpenSees.VectorWrapper(new double[]{ coreY, coreZ }), new OpenSees.VectorWrapper(new double[]{ coreY, ncoreZ})),

				//# Define remaining steel in the y direction
				// layer straight $steelID $numBars $barArea [expr $coreY-$spacingY] $coreZ [expr $ncoreY+$spacingY] $coreZ
				//new OpenSees.Materials.Sections.Repres.StraightReinfLayerWrapper(rebarMat.Key,numYBars,barArea,new OpenSees.VectorWrapper(new double[]{ coreY -spacingY, coreZ }), new OpenSees.VectorWrapper(new double[]{ ncoreY +spacingY, coreZ})),
				//layer straight $steelID $numBars $barArea [expr $coreY-$spacingY] $ncoreZ [expr $ncoreY+$spacingY] $ncoreZ
				//new OpenSees.Materials.Sections.Repres.StraightReinfLayerWrapper(rebarMat.Key,numYBars,barArea,new OpenSees.VectorWrapper(new double[]{ coreY -spacingY, ncoreZ }), new OpenSees.VectorWrapper(new double[]{ ncoreY +spacingY, ncoreZ})),
	   }
	);

	var matDic = new Dictionary<int, OpenSees.Materials.Uniaxials.UniaxialMaterialWrapper>
	{
		{ rebarMat.Key, rebarMat.Value },
		{ coverMat.Key, coverMat.Value },
		{ coreMat.Key, coreMat.Value }
	};
	var theSection = new OpenSees.Materials.Sections.FiberSection3dWrapper(1, repres, matDic,
		new Materials.Sections.Repres.UniaxialFiber3dWrapper[] {
			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ coreY, coreZ})),
			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ 0, coreZ})),
			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ -coreY, coreZ})),

			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ coreY, 0})),
			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ -coreY, 0})),

			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ coreY, -coreZ})),
			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ 0, -coreZ})),
			new Materials.Sections.Repres.UniaxialFiber3dWrapper(1,rebarMat.Value,barArea,new VectorWrapper(new double[]{ -coreY, -coreZ})),

		}, GJ);
	return theSection;
}
                        

wipe
# Units: kip, in, sec

# ----------------------------
# Start of model generation
# ----------------------------

# Create ModelBuilder with 3 dimensions and 6 DOF/node
model BasicBuilder -ndm 3 -ndf 6

# Define geometry
# ---------------

# Set parameters for model geometry
set h 144.0;       # Story height
set by 240.0;      # Bay width in Y-direction
set bx 240.0;      # Bay width in X-direction

# Create nodes
#    tag             X             Y           Z 
node  1  [expr -$bx/2] [expr  $by/2]           0
node  2  [expr  $bx/2] [expr  $by/2]           0
node  3  [expr  $bx/2] [expr -$by/2]           0 
node  4  [expr -$bx/2] [expr -$by/2]           0 

node  5  [expr -$bx/2] [expr  $by/2]          $h 
node  6  [expr  $bx/2] [expr  $by/2]          $h 
node  7  [expr  $bx/2] [expr -$by/2]          $h 
node  8  [expr -$bx/2] [expr -$by/2]          $h 

node 10  [expr -$bx/2] [expr  $by/2] [expr 2*$h]
node 11  [expr  $bx/2] [expr  $by/2] [expr 2*$h] 
node 12  [expr  $bx/2] [expr -$by/2] [expr 2*$h] 
node 13  [expr -$bx/2] [expr -$by/2] [expr 2*$h] 

node 15  [expr -$bx/2] [expr  $by/2] [expr 3*$h] 
node 16  [expr  $bx/2] [expr  $by/2] [expr 3*$h] 
node 17  [expr  $bx/2] [expr -$by/2] [expr 3*$h] 
node 18  [expr -$bx/2] [expr -$by/2] [expr 3*$h]

# Master nodes for rigid diaphragm
#    tag X Y          Z 
node  9  0 0         $h 
node 14  0 0 [expr 2*$h]
node 19  0 0 [expr 3*$h]

# Set base constraints
#   tag DX DY DZ RX RY RZ
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

# Define rigid diaphragm multi-point constraints
#               normalDir  master     slaves
rigidDiaphragm     3          9     5  6  7  8
rigidDiaphragm     3         14    10 11 12 13
rigidDiaphragm     3         19    15 16 17 18

# Constraints for rigid diaphragm master nodes
#   tag DX DY DZ RX RY RZ
fix  9   0  0  1  1  1  0
fix 14   0  0  1  1  1  0
fix 19   0  0  1  1  1  0

# Define materials for nonlinear columns
# --------------------------------------
# CONCRETE
# Core concrete (confined)
#                           tag  f'c  epsc0  f'cu  epscu
uniaxialMaterial Concrete01  1  -5.0 -0.005  -3.5  -0.02

# Cover concrete (unconfined)
set fc 4.0
uniaxialMaterial Concrete01  2   -$fc -0.002   0.0 -0.006

# STEEL
# Reinforcing steel
#                        tag fy   E     b
uniaxialMaterial Steel01  3  60 30000 0.02

# Column width
set h 18.0

# Source in a procedure for generating an RC fiber section
source RCsection.tcl

# Column torsional stiffness
set GJ 1.0e10;

# Call the procedure to generate the column section
#          id  h  b cover core cover steel nBars barArea nfCoreY nfCoreZ nfCoverY nfCoverZ GJ
RCsection   1 $h $h   2.5    1    2     3     3    0.79       8       8       10       10 $GJ

# Concrete elastic stiffness
set E [expr 57000.0*sqrt($fc*1000)/1000];



# Linear elastic torsion for the column
# uniaxialMaterial Elastic 10 $GJ

# Attach torsion to the RC column section
#                 tag uniTag uniCode       secTag
# section Aggregator 2    10      T    -section 1


set colSec 1

# Define column elements
# ----------------------

#set PDelta "ON"
set PDelta "OFF"

# Geometric transformation for columns
if {$PDelta == "ON"} {
   #                           tag  vecxz
   geomTransf LinearWithPDelta  1   1 0 0
} else {
   geomTransf Linear  1   1 0 0
}

# Number of column integration points (sections)


set np 4

# Create the nonlinear column elements
#                           tag ndI ndJ nPts   secID  transf
element nonlinearBeamColumn  1   1   5   $np  $colSec    1
element nonlinearBeamColumn  2   2   6   $np  $colSec    1
element nonlinearBeamColumn  3   3   7   $np  $colSec    1
element nonlinearBeamColumn  4   4   8   $np  $colSec    1

element nonlinearBeamColumn  5   5  10   $np  $colSec    1
element nonlinearBeamColumn  6   6  11   $np  $colSec    1
element nonlinearBeamColumn  7   7  12   $np  $colSec    1
element nonlinearBeamColumn  8   8  13   $np  $colSec    1

element nonlinearBeamColumn  9  10  15   $np  $colSec    1
element nonlinearBeamColumn 10  11  16   $np  $colSec    1
element nonlinearBeamColumn 11  12  17   $np  $colSec    1
element nonlinearBeamColumn 12  13  18   $np  $colSec    1

# Define beam elements
# --------------------

# Define material properties for elastic beams
# Using beam depth of 24 and width of 18
# --------------------------------------------
set Abeam [expr 18*24];
# "Cracked" second moments of area
set Ibeamzz [expr 0.5*1.0/12*18*pow(24,3)];
set Ibeamyy [expr 0.5*1.0/12*24*pow(18,3)];

# Define elastic section for beams
#               tag  E    A      Iz       Iy     G   J
section Elastic  3  $E $Abeam $Ibeamzz $Ibeamyy $GJ 1.0
set beamSec 3

# Geometric transformation for beams
#                tag  vecxz
geomTransf Linear 2   1 1 0

# Number of beam integration points (sections)
set np 3

# Create the beam elements
#                           tag ndI ndJ nPts    secID   transf
element nonlinearBeamColumn  13  5   6   $np  $beamSec     2
element nonlinearBeamColumn  14  6   7   $np  $beamSec     2
element nonlinearBeamColumn  15  7   8   $np  $beamSec     2
element nonlinearBeamColumn  16  8   5   $np  $beamSec     2

element nonlinearBeamColumn  17 10  11   $np  $beamSec     2
element nonlinearBeamColumn  18 11  12   $np  $beamSec     2
element nonlinearBeamColumn  19 12  13   $np  $beamSec     2
element nonlinearBeamColumn  20 13  10   $np  $beamSec     2

element nonlinearBeamColumn  21 15  16   $np  $beamSec     2
element nonlinearBeamColumn  22 16  17   $np  $beamSec     2
element nonlinearBeamColumn  23 17  18   $np  $beamSec     2
element nonlinearBeamColumn  24 18  15   $np  $beamSec     2

# Define gravity loads
# --------------------
# Gravity load applied at each corner node
# 10% of column capacity
set p [expr 0.1*$fc*$h*$h]

# Mass lumped at master nodes
set g 386.4;            # Gravitational constant
set m [expr (4*$p)/$g]

# Rotary inertia of floor about master node
set i [expr $m*($bx*$bx+$by*$by)/12.0]
#puts "mass: $m $m  0  0  0 $i"
# Set mass at the master nodes
#    tag MX MY MZ RX RY RZ
mass  9  $m $m  0  0  0 $i
mass 14  $m $m  0  0  0 $i
mass 19  $m $m  0  0  0 $i

# Define gravity loads
pattern Plain 1 Constant {
   foreach node {5 6 7 8  10 11 12 13  15 16 17 18} {
      load $node 0.0 0.0 -$p 0.0 0.0 0.0
   }
}

#set rayleigh damping factors
rayleigh 0.0 0.0 0.0 0.0018


# Define earthquake excitation
# ----------------------------
# Set up the acceleration records for Tabas fault normal and fault parallel
set tabasFN "Path -filePath tabasFN.txt -dt 0.02 -factor $g"
set tabasFP "Path -filePath tabasFP.txt -dt 0.02 -factor $g"

# Define the excitation using the Tabas ground motion records
#                         tag dir         accel series args
pattern UniformExcitation  2   1  -accel      $tabasFN
pattern UniformExcitation  3   2  -accel      $tabasFP


# ----------------------- 
# End of model generation
# -----------------------


# ----------------------------
# Start of analysis generation
# ----------------------------

# Create the convergence test
#                tol   maxIter  printFlag
test EnergyIncr 1.0e-8   20         3

# Create the solution algorithm
algorithm Newton

# Create the constraint handler
constraints Transformation

# Create the time integration scheme
#                   gamma beta
integrator Newmark   0.5  0.25 

# Create the system of equation storage and solver
system SparseGeneral -piv

# Create the DOF numberer
numberer Plain

# Create the transient analysis
analysis Transient

# --------------------------
# End of analysis generation
# --------------------------
# ----------------------------

# Record DOF 1 and 2 displacements at nodes 9, 14, and 19
recorder Node -file Node51.out -time -node 9 14 19 -dof 1 2 disp
recorder plot Node51.out Node9_14_19_Xdisp 10 340 300 300 -columns 1 2 -columns 1 4 -columns 1 6  -dT 1.0

# --------------------------
# End of recorder generation
# --------------------------


# --------------------
# Perform the analysis
# --------------------

# Analysis duration of 20 seconds
#       numSteps  dt

puts "[eigen -fullGenLapack 3]";
set ok [analyze   2000   0.01]
if {$ok != 0} {
    puts "analysis FAILED"
} else {
    puts "analysis SUCCESSFULL"
}
#http://opensees.berkeley.edu/OpenSees/manuals/ExamplesManual/HTML/871.htm

# Define a procedure which generates a rectangular reinforced concrete section
# with one layer of steel evenly distributed around the perimeter and a confined core.
# 
#                       y
#                       |
#                       |
#                       |    
#             ---------------------
#             |\                 /|
#             | \---------------/ |
#             | |               | |
#             | |               | |
#  z ---------| |               | |  h
#             | |               | |
#             | |               | |
#             | /---------------\ |
#             |/                 \|
#             ---------------------
#                       b
#
# Formal arguments
#    id - tag for the section that is generated by this procedure
#    h - overall height of the section (see above)
#    b - overall width of the section (see above)
#    cover - thickness of the cover patches
#    coreID - material tag for the core patch
#    coverID - material tag for the cover patches
#    steelID - material tag for the reinforcing steel
#    numBars - number of reinforcing bars on any given side of the section
#    barArea - cross-sectional area of each reinforcing bar
#    nfCoreY - number of fibers in the core patch in the y direction
#    nfCoreZ - number of fibers in the core patch in the z direction
#    nfCoverY - number of fibers in the cover patches with long sides in the y direction
#    nfCoverZ - number of fibers in the cover patches with long sides in the z direction
#
# Notes
#    The thickness of cover concrete is constant on all sides of the core.
#    The number of bars is the same on any given side of the section.
#    The reinforcing bars are all the same size.
#    The number of fibers in the short direction of the cover patches is set to 1.
# 
proc RCsection {id h b cover coreID coverID steelID numBars barArea nfCoreY nfCoreZ nfCoverY nfCoverZ GJ} {

   # The distance from the section z-axis to the edge of the cover concrete
   # in the positive y direction
   set coverY [expr $h/2.0]

   # The distance from the section y-axis to the edge of the cover concrete
   # in the positive z direction
   set coverZ [expr $b/2.0]

   # The negative values of the two above
   set ncoverY [expr -$coverY]
   set ncoverZ [expr -$coverZ]

   # Determine the corresponding values from the respective axes to the
   # edge of the core concrete
   set coreY [expr $coverY-$cover]
   set coreZ [expr $coverZ-$cover]
   set ncoreY [expr -$coreY]
   set ncoreZ [expr -$coreZ]

   # Define the fiber section
   section fiberSec $id -GJ $GJ {

	# Define the core patch
	patch quadr $coreID $nfCoreZ $nfCoreY $ncoreY $coreZ $ncoreY $ncoreZ $coreY $ncoreZ $coreY $coreZ
      
	# Define the four cover patches
	patch quadr $coverID 1 $nfCoverY $ncoverY $coverZ $ncoreY $coreZ $coreY $coreZ $coverY $coverZ
	patch quadr $coverID 1 $nfCoverY $ncoreY $ncoreZ $ncoverY $ncoverZ $coverY $ncoverZ $coreY $ncoreZ
	patch quadr $coverID $nfCoverZ 1 $ncoverY $coverZ $ncoverY $ncoverZ $ncoreY $ncoreZ $ncoreY $coreZ
	patch quadr $coverID $nfCoverZ 1 $coreY $coreZ $coreY $ncoreZ $coverY $ncoverZ $coverY $coverZ

	#fiber  $ncoreY $coreZ $barArea $steelID 
	#fiber  0 $coreZ $barArea $steelID 
	#fiber  $coreY $coreZ $barArea $steelID 

	#fiber  $ncoreY 0 $barArea $steelID 
	#fiber  $coreY 0 $barArea $steelID 

	#fiber  $ncoreY $ncoreZ $barArea $steelID 
	#fiber  0 $ncoreZ $barArea $steelID 
	#fiber  $coreY $ncoreZ $barArea $steelID 

	# Define the steel along constant values of y (in the z direction)
	layer straight $steelID $numBars $barArea $ncoreY $coreZ $ncoreY $ncoreZ
	layer straight $steelID $numBars $barArea $coreY $coreZ $coreY $ncoreZ

	# Determine the spacing for the remaining bars in the y direction
	set spacingY [expr ($coreY-$ncoreY)/($numBars-1)]

	# Avoid double counting bars
	set numBars [expr $numBars-2]
	# Define remaining steel in the y direction
	layer straight $steelID $numBars $barArea [expr $coreY-$spacingY] $coreZ [expr $ncoreY+$spacingY] $coreZ
	layer straight $steelID $numBars $barArea [expr $coreY-$spacingY] $ncoreZ [expr $ncoreY+$spacingY] $ncoreZ
   }

}
                            

Reference : http://opensees.berkeley.edu/OpenSees/manuals/ExamplesManual/HTML/871.htm

Ground-motion file  : tabasFN.txt

Ground-motion file  : tabasFP.txt