# Input file for Analysis of Highway Bridge Structure.

You will need Aladdin 2.0 to run this problem.

```/*
*  =====================================================================
*  Analysis of 4 span isolated bridge, 3 piers, 2 abutments, 5 isolators
*
*  Modeling Assumptions:
*
*     --  Bottoms of the abuntments' isolators and piers are fixed.
*     --  Deck is modeled as FRAME_3D elements.
*     --  Pier columns are modeled as FRAME_3D elements.
*     --  Isolators are modeled as FIBER_3D elements.
*     --  Perfect elastic-plastic for isolators.
*     --  Peak earthquake = 0.5*g. (rec.6)
*     --  5% of structural damping is added, C=A0*M+A1+K.
*
*  Written By : Wane-Jang Lin                               October 1997
*  =====================================================================
*/

NDimension         = 3;
NDofPerNode        = 6;
MaxNodesPerElement = 2;
GaussIntegPts      = 2;

StartMesh();

H2    = 14 m;
H3    =  7 m;
H4    = 21 m;
iso_h = 20 cm;

/* Assemble mesh of nodes and elements */

total_node = 0;
x = 0 m;  y = 0 m;  z = 0 m;
for( node_no=total_node+1 ; x<=200m ; node_no=node_no+1 ) {
x = x + 50 m;
}

total_node = node_no-1;
x = 0 m;  y = -iso_h;  z = 0 m;
for( node_no=total_node+1 ; x<=200m ; node_no=node_no+1 ) {
x = x + 50 m;
}

total_node = node_no-1;
x = 50 m;  y = -iso_h-3.5m;  z = 0 m;
for( node_no=total_node+1 ; y>=-(iso_h+H2) ; node_no=node_no+1 ) {
y = y - 3.5 m;
}

total_node = node_no-1;
fix2       = total_node;
x = 100 m;  y = -iso_h-3.5m;  z = 0 m;
for( node_no=total_node+1 ; y>=-(iso_h+H3) ; node_no=node_no+1 ) {
y = y - 3.5 m;
}

total_node = node_no-1;
fix3       = total_node;
x = 150 m;  y = -iso_h-3.5m;  z = 0 m;
for( node_no=total_node+1 ; y>=-(iso_h+H4) ; node_no=node_no+1 ) {
y = y - 3.5 m;
}

total_node = node_no-1;
fix4       = total_node;

total_elmt = 0;
for( elmt_no=total_elmt+1 ; elmt_no<=4 ; elmt_no=elmt_no+1 ) {
AddElmt( elmt_no, [elmt_no, elmt_no+1], "deck_attr" );
}

total_elmt = 9;
node_no = fix2 - H2/(3.5m) + 1;
AddElmt( total_elmt+1, [7, node_no], "column_attr2" );

total_elmt = total_elmt+1;
for( elmt_no=total_elmt+1 ; node_no < fix2 ; elmt_no=elmt_no+1 ) {
AddElmt( elmt_no, [node_no, node_no+1], "column_attr2" );
node_no = node_no+1;
}

total_elmt = elmt_no-1;
bot_elmt2  = total_elmt;
node_no    = fix3 - H3/(3.5m) + 1;
AddElmt( total_elmt+1, [8, node_no], "column_attr3" );
total_elmt = total_elmt+1;
for( elmt_no=total_elmt+1 ; node_no < fix3 ; elmt_no=elmt_no+1 ) {
AddElmt( elmt_no, [node_no, node_no+1], "column_attr3" );
node_no = node_no+1;
}

total_elmt = elmt_no-1;
bot_elmt3  = total_elmt;
node_no    = fix4 - H4/(3.5m) + 1;
AddElmt( total_elmt+1, [9, node_no], "column_attr4" );
total_elmt = total_elmt+1;
for( elmt_no=total_elmt+1 ; node_no < fix4 ; elmt_no=elmt_no+1 ) {
AddElmt( elmt_no, [node_no, node_no+1], "column_attr4" );
node_no = node_no+1;
}

total_elmt = elmt_no-1;
bot_elmt4  = total_elmt;

ElementAttr("deck_attr") { type     = "FRAME_3D";
section  = "deck_sec";
material = "deck_mat";
}
SectionAttr("deck_sec")  { Iyy = 103.77 m^4; Izz = 5.26 m^4;
J       = 16.13 m^4;
area    = 6.88 m^2;
unit_weight = 200 kN/m;
}
MaterialAttr("deck_mat") { poisson = 0.25;
E       = 2.3e7 kN/m^2;
}

b  = 2.2    m;  h = 4.0    m;
Ix = 7.39 m^4; Iy = 0.67 m^4;
K2 =       30700 kN/m; K3 =      124400 kN/m; K4 = 13650 kN/m;
E2 = K2*H2*H2*H2/3/Ix; E3 = K3*H3*H3*H3/3/Ix; E4 = K4*H4*H4*H4/3/Ix;

ElementAttr("column_attr2") { type     = "FRAME_3D";
section  = "col_sec";
material = "col_mat2";
}

ElementAttr("column_attr3") { type     = "FRAME_3D";
section  = "col_sec";
material = "col_mat3";
}

ElementAttr("column_attr4") { type     = "FRAME_3D";
section  = "col_sec";
material = "col_mat4";
}

SectionAttr("col_sec") { area  = 4.16 m^2;
width = h;  depth = b;
Iyy   = Ix; Izz   = Iy;
unit_weight  = 104 kN/m;
shear_factor = 1.0;
}

MaterialAttr("col_mat2") { poisson = 0.25;  E = E2; }
MaterialAttr("col_mat3") { poisson = 0.25;  E = E3; }
MaterialAttr("col_mat4") { poisson = 0.25;  E = E4; }

h = 50 cm;
b = 50 cm;
G1  = 4240 kN/m^2;  G2 = 12880 kN/m^2;   G3 = 9200 kN/m^2;   G4 = 37600 kN/m^2;
fvy = 8400 kN/m^2;

ElementAttr("isolator_attr1") { type     = "FIBER_3D";
section  = "iso_sec";
material = "iso_mat1";
fiber    = "iso_fib";
}

ElementAttr("isolator_attr2") { type     = "FIBER_3D";
section  = "iso_sec";
material = "iso_mat2";
fiber    = "iso_fib";
}

ElementAttr("isolator_attr3") { type     = "FIBER_3D";
section  = "iso_sec";
material = "iso_mat3";
fiber    = "iso_fib";
}

ElementAttr("isolator_attr4") { type     = "FIBER_3D";
section  = "iso_sec";
material = "iso_mat4";
fiber    = "iso_fib";
}

SectionAttr("iso_sec") { area  = b*h;
width = b;
depth = h;
unit_weight  = 0.0001 kN/m;
shear_factor = 1.0;
}

MaterialAttr("iso_mat1") { poisson = 0.25; G = G1; Gt = 0.001*G1; shear_yield = fvy/2; }
MaterialAttr("iso_mat2") { poisson = 0.25; G = G2; Gt = 0.001*G2; shear_yield = fvy; }
MaterialAttr("iso_mat3") { poisson = 0.25; G = G3; Gt = 0.001*G3; shear_yield = fvy; }
MaterialAttr("iso_mat4") { poisson = 0.25; G = G4; Gt = 0.001*G4; shear_yield = fvy; }

no_fiber = 4;
iso_fcoord = Matrix([ 2, no_fiber ]);
iso_farea  = Matrix([ 1, no_fiber ]);
iso_fmap   = Matrix([ 1, no_fiber ]);

iso_fcoord[1][1] =  b/4;
iso_fcoord[2][1] =  h/4;
iso_fcoord[1][2] = -b/4;
iso_fcoord[2][2] =  h/4;
iso_fcoord[1][3] =  b/4;
iso_fcoord[2][3] = -h/4;
iso_fcoord[1][4] = -b/4;
iso_fcoord[2][4] = -h/4;

for( i=1 ; i <= no_fiber ; i=i+1 ) {
iso_farea[1][i]  = b*h/no_fiber;
iso_fmap[1][i]   = 1;
}

iso_fattr = Matrix([ 3, 1 ]);
iso_fattr[3][1] = fvy*1e6;

FiberAttr( no_fiber, "iso_fib" ) { FiberMaterialAttr = iso_fattr;
FiberCoordinate   = iso_fcoord;
FiberArea         = iso_farea;
FiberMaterialMap  = iso_fmap;
}

/* Apply boundary conditions */

FixNode(   1, [1,0,0,1,0,0]);
FixNode(   5, [1,0,0,1,0,0]);
FixNode(   6, [1,1,1,1,1,1]);
FixNode(  10, [1,1,1,1,1,1]);
FixNode(fix2, [1,1,1,1,1,1]);
FixNode(fix3, [1,1,1,1,1,1]);
FixNode(fix4, [1,1,1,1,1,1]);

/* Compile and print f.e. mesh    */

EndMesh();
PrintMesh();

/* Compute initial stiffness and mass matrices */

stiff = Stiff();
mass  = Mass([1]);
mass_inv = Inverse( mass );

/* Compute and print eigenvalues */

no_eigen   = 2;
eigen      = Eigen(stiff, mass, [no_eigen]);
eigenvalue = Eigenvalue(eigen);
for(i=1 ; i <= no_eigen ; i=i+1) {
w = sqrt( eigenvalue[i][1] );
print "Mode",i,",\tw=",w,",\tperiod=",2*PI/w,"\n";
}

eigenvector = Eigenvector(eigen);
PrintMatrix(eigenvector);

wa = sqrt( eigenvalue[1][1] );
wb = sqrt( eigenvalue[2][1] );

/* Setup Rayleigh damping and damping matrix */

rdamping = 0.20;
A0   = 2*rdamping*wa*wb/(wa+wb);
A1   = 2*rdamping/(wa+wb);
damp = A0*mass + A1*stiff;

/* El Centro Earthquake record */

Elcentro1 = ColumnUnits( [
4.56;    13.77;     6.13;     3.73;     1.32;    -6.81;   -16.22;   -22.41;
-26.67;   -24.10;   -21.86;   -14.02;    -1.62;     1.67;    18.12;    39.26;
28.40;    16.72;    16.06;    12.54;     0.15;    -8.60;   -14.00;   -12.01;
-3.34;     5.14;    14.04;    16.36;    29.90;    43.09;    32.13;    17.28;
18.44;    22.31;    13.24;    -0.91;    -6.81;   -15.20;   -14.88;    -9.54;
-13.32;   -10.60;    11.80;    10.20;   -12.90;   -22.88;     1.42;    21.73;
7.92;     6.47;    11.39;    -8.29;   -33.44;   -26.25;   -23.99;   -31.05;
-25.42;   -26.97;   -38.48;   -35.01;   -27.78;   -33.35;   -29.56;   -13.60;
-13.97;   -11.63;    -1.31;    13.86;    35.71;    27.62;    25.33;    46.45;
44.12;    28.02;    19.13;    12.57;    11.60;    12.29;     9.73;    10.38;
18.93;    29.43;    30.25;    26.23;    23.15;    22.29;    21.11;    14.84;
11.02;    14.74;    21.65;    26.66;    28.45;    30.34;    30.21;    28.82;
31.31;    26.66;    15.29;     9.10;     3.58;    -6.65;   -15.33;   -18.05;
-17.46;   -11.68;    -3.48;    -5.25;   -15.43;   -25.38;   -32.32;   -29.67;
-22.01;   -17.97;   -22.34;   -30.51;   -38.89;   -45.29;   -55.26;   -67.84;
-76.74;   -84.58;   -86.83;   -86.63;   -80.87;   -69.83;   -59.81;   -54.51;
-48.94;   -39.94;   -38.67;   -40.13;   -38.04;   -34.78;   -31.30;   -21.29;
-8.85;     3.58;    13.97;    30.88;    35.09;     8.84;   -10.41;   -16.30;
-20.10;    -3.69;    33.71;    51.69;    43.52;    27.04;    34.70;    38.30;
36.60;    29.85;    -9.22;   -23.98;   -17.28;   -14.20;    11.58;    47.99;
64.26;    70.36;    68.93;    62.67;    45.30;    28.43;    23.73;    28.51;
34.46;    43.12;    43.39;    26.64;    21.23;    31.23;    40.30;    16.76;
-10.80;   -22.88;    -0.51;    12.51;   -10.01;   -12.52;    -6.05;   -24.18;
-43.17;     7.01;    50.06;    54.19;    68.83;    69.26;    57.53;     8.73;
-38.98;   -43.59;   -51.12;   -41.86;   -19.51;   -15.78;   -11.54;    -2.18;
-2.19;    20.84;    23.40;    -0.91;   -11.45;   -22.94;   -40.84;   -53.17;
-44.02;   -21.96;    -3.06;    -0.49;    -9.59;   -11.83;    -9.86;    -7.02;
0.88;    24.00;    50.57;    56.01;    51.75;    52.30;    47.79;    16.77;
-6.09;   -20.13;   -20.42;    -4.35;    -0.33;    -6.33;     1.74;    19.90;
37.97;    53.36;    57.68;    53.32;    28.37;     7.54;     6.90;     4.73;
-4.97;   -21.97;   -28.16;   -23.25;   -28.43;   -18.57;    -3.93;   -16.60;
-31.92;   -35.08;   -34.77;   -30.85;   -30.37;    -6.83;    24.57;    27.39;
23.72;    29.30;     5.79;   -35.35;   -54.31;   -57.00;   -59.38;   -19.00;
32.14;     7.01;   -21.44;   -21.26;   -26.01;   -35.66;   -33.68;   -24.09;
-21.39;   -26.29;   -27.73;   -23.91;   -14.27;    -6.41;    -0.59;     4.21;
1.06;    -4.59;   -15.16;   -18.12;    -7.09;     2.41;     6.08;     1.20;
-10.53;    -8.69;    -4.70;   -13.05;    -9.51;    -0.67;    -7.10;    -7.19;
5.36;     9.80;     7.87;    12.63;    20.28;    22.44;    24.86;    24.31;
15.25;     6.54;     6.86;    11.08;    18.12;    23.64;    22.47;    23.76;
21.57;    19.80;    25.47;    27.94;    14.81;    -6.94;   -13.57;    -7.80;
-2.65;    -3.97;    -1.01;     0.20;     5.14;     2.31;   -16.89;   -18.18;
-11.55;   -12.10;    -3.07;     6.16;     8.87;    19.72;    21.31;    16.45;
12.70;     4.15;    -1.59;     0.57;    11.04;    21.39;    27.42;    15.82;
-9.25;   -23.67;   -17.34;     0.67;     8.69;    11.01;    12.01;    11.12;
13.62;    11.31;     9.71;    17.68;    15.83;     5.59;     0.60;    -2.53;
-4.46;    -4.09;    -7.81;   -15.06;   -19.43;   -15.49;    -8.93;    -5.44;
-1.33;    -1.35;     7.19;    12.52;     0.21;    -4.89;     0.22;    11.51;
20.17;    18.62;    16.72;    13.74;    10.62;    11.14;    12.74;    13.34;
14.16;    14.59;    13.81;    10.94;     5.91;     5.74;    14.53;    19.56;
9.27;    -6.75;   -14.69;   -14.25;   -11.48;    -3.96;     6.54;    10.59;
2.26;   -12.97;   -24.01;   -25.73;   -18.75;    -6.99;    -2.22;     0.17;
3.63;    -1.70;    -4.73;     6.13;    20.26;    22.36;    19.58;     9.14;
-0.09;    -0.23;    -1.00;     4.82;    15.39;    19.75;     7.45;    -0.99;
0.77;    -1.64;    -9.34;   -14.34;   -13.86;    -6.61;     5.49;     6.30;
2.08;    -2.88;    -2.25;     2.99;     3.30;    -1.57;    -6.86;    -8.44;
-11.03;   -16.68;   -17.43;   -12.39;   -12.20;   -15.43;   -16.28;   -17.99;
-20.21;   -21.58;   -23.05;   -23.79;   -23.64;   -21.60;   -15.32;    -5.72;
11.14;    28.30;    33.22;    30.27;    25.49;    18.45;     7.62;    -1.66;
-7.86;   -10.04;    -5.71;    -5.16;   -12.56;   -19.39;   -22.16;   -18.18;
-12.00;    -3.60;     3.99;     3.89;     0.85;     4.01;     6.05;     0.68;
-0.14;     1.50;     3.17;     4.96;    -0.03;    -7.20;    -4.85;    -0.60;
-3.05;    -6.01;    -6.08;    -6.40;    -6.70;    -5.70;    -6.14;    -8.41;
-10.05;   -12.35;   -15.72;     0.00
], [cm/sec/sec] );

dimen = Dimension(Elcentro1);
x     = dimen[1][1];

/* Interpolate accelerations */

divid_no = 16;
Elcentro = ColumnUnits( Matrix([ divid_no*x, 1 ]),[in/sec/sec] );
Elcentro[1][1] = Elcentro1[1][1]/16;
Elcentro[2][1] = Elcentro1[1][1]/8;
Elcentro[3][1] = Elcentro1[1][1]*3/16;
Elcentro[4][1] = Elcentro1[1][1]/4;
Elcentro[5][1] = Elcentro1[1][1]*5/16;
Elcentro[6][1] = Elcentro1[1][1]*3/8;
Elcentro[7][1] = Elcentro1[1][1]*7/16;
Elcentro[8][1] = Elcentro1[1][1]/2;
Elcentro[9][1] = Elcentro1[1][1]*9/16;
Elcentro[10][1]= Elcentro1[1][1]*5/8;
Elcentro[11][1]= Elcentro1[1][1]*11/16;
Elcentro[12][1]= Elcentro1[1][1]*3/4;
Elcentro[13][1]= Elcentro1[1][1]*13/16;
Elcentro[14][1]= Elcentro1[1][1]*7/8;
Elcentro[15][1]= Elcentro1[1][1]*15/16;
Elcentro[divid_no*x][1] = Elcentro1[x][1];
for( i=1 ; i < x; i=i+1 ) {
Elcentro[divid_no*i][1] = Elcentro1[i][1];
Elcentro[divid_no*i+1][1] = ( Elcentro1[i][1]*15+Elcentro1[i+1][1] )/16;
Elcentro[divid_no*i+2][1] = ( Elcentro1[i][1]*7 +Elcentro1[i+1][1] )/8;
Elcentro[divid_no*i+3][1] = ( Elcentro1[i][1]*13+Elcentro1[i+1][1]*3 )/16;
Elcentro[divid_no*i+4][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1] )/4;
Elcentro[divid_no*i+5][1] = ( Elcentro1[i][1]*11+Elcentro1[i+1][1]*5 )/16;
Elcentro[divid_no*i+6][1] = ( Elcentro1[i][1]*5 +Elcentro1[i+1][1]*3 )/8;
Elcentro[divid_no*i+7][1] = ( Elcentro1[i][1]*9 +Elcentro1[i+1][1]*7 )/16;
Elcentro[divid_no*i+8][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1] )/2;
Elcentro[divid_no*i+9][1] = ( Elcentro1[i][1]*7 +Elcentro1[i+1][1]*9 )/16;
Elcentro[divid_no*i+10][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1]*5 )/8;
Elcentro[divid_no*i+11][1] = ( Elcentro1[i][1]*5 +Elcentro1[i+1][1]*11 )/16;
Elcentro[divid_no*i+12][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1]*3 )/4;
Elcentro[divid_no*i+13][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1]*13 )/16;
Elcentro[divid_no*i+14][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1]*7 )/8;
Elcentro[divid_no*i+15][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1]*15 )/16;
}

/*
divid_no = 8;
Elcentro = ColumnUnits( Matrix([ divid_no*x, 1 ]),[in/sec/sec] );
Elcentro[1][1] = Elcentro1[1][1]/8;
Elcentro[2][1] = Elcentro1[1][1]/4;
Elcentro[3][1] = Elcentro1[1][1]*3/8;
Elcentro[4][1] = Elcentro1[1][1]/2;
Elcentro[5][1] = Elcentro1[1][1]*5/8;
Elcentro[6][1] = Elcentro1[1][1]*3/4;
Elcentro[7][1] = Elcentro1[1][1]*7/8;
Elcentro[divid_no*x][1] = Elcentro1[x][1];
for( i=1 ; i < x; i=i+1 ) {
Elcentro[divid_no*i][1] = Elcentro1[i][1];
Elcentro[divid_no*i+1][1] = ( Elcentro1[i][1]*7+Elcentro1[i+1][1] )/8;
Elcentro[divid_no*i+2][1] = ( Elcentro1[i][1]*3+Elcentro1[i+1][1] )/4;
Elcentro[divid_no*i+3][1] = ( Elcentro1[i][1]*5+Elcentro1[i+1][1]*3 )/8;
Elcentro[divid_no*i+4][1] = ( Elcentro1[i][1]  +Elcentro1[i+1][1] )/2;
Elcentro[divid_no*i+5][1] = ( Elcentro1[i][1]*3+Elcentro1[i+1][1]*5 )/8;
Elcentro[divid_no*i+6][1] = ( Elcentro1[i][1]  +Elcentro1[i+1][1]*3 )/4;
Elcentro[divid_no*i+7][1] = ( Elcentro1[i][1]  +Elcentro1[i+1][1]*7 )/8;
}
*/

x = Max( [ Max(Elcentro), abs(Min(Elcentro)) ] );
Elcentro = (0.5*9.81/x)*Elcentro;

print "\n*** Time History Analysis :";
print "\n*** Use Average Acceleration Method :\n\n";

/* ----------------------------------------------------- */
/* Setup Time History Analysis  (and initial conditions) */
/* ----------------------------------------------------- */

NodeLoad( 1, [ 0 kN, 0 kN, 0 kN, 0 kN*m, 0 kN*m, 0 kN*m] );
P_old = P_ext;

/* Setup initial displacement, velocity and acceleration to zeros */

displ    = Solve( stiff, P_ext );
velocity = displ/(1 sec);
accel    = velocity/(1 sec);

new_displ    = displ;
new_velocity = velocity;
new_accel    = accel;
displ_i      = displ;

/* Setup initial internal force and damping force */

Fd_old = damp*velocity;
Fs_i   = Fs_old;

/* Setup the influence vector */

r = displ/(1 m);
accel_dir = 3;
for( i=1 ; i<=total_node ; i=i+1 ) {
dof = GetDof([i]);
if( dof[1][accel_dir] > 0 ) {
r[ dof[1][accel_dir] ][1] = 1;
}
}

/* Setup time interval and analysis time */

dt     = (0.02 sec)/divid_no;
time   = 0.0  sec;
tol    = 0.0001;
stepno = 0;

dimen = Dimension(Elcentro);
total_stepno = dimen[1][1] + 250*divid_no;
total_time   = total_stepno * dt;
quake_time   = dimen[1][1] * dt;
print "        dt =",dt,"\n";
print "total step =",total_stepno," time       =",total_time,"\n";
print "quake step =", dimen[1][1]," quake time =",quake_time,"\n";

/* Allocate response matrices */

deck_displ     = ColumnUnits(Zero([total_stepno,5]),[m]);
pier_displ     = ColumnUnits(Zero([total_stepno,3]),[m]);
bottom_shear   = ColumnUnits(Zero([total_stepno,3]),[N]);
isolator_shear = ColumnUnits(Zero([total_stepno,5]),[N]);

energy  = ColumnUnits(Zero([total_stepno,4]),[Jou]);

/* Time-History Analysis */

while(stepno < total_stepno) {

/* Update time and step no */

time   = time + dt;
stepno = stepno + 1;

print "*** Start at step  ", stepno, " : TIME = ", time, "\n";

if( stepno <= dimen[1][1] ) then {
P_ext = -mass*r*Elcentro[stepno][1];
} else {
P_ext = -mass*r*(0.0 m/sec/sec);
}

dPeff = P_ext - P_old + ((4/dt)*mass + 2*damp)*velocity + 2*mass*accel;

/* while loop to check converge, Keff*U = Peff */

dp  = displ - displ;
err = 1 + tol;

ii = 1;
while( err > tol ) {

/* Compute effective stiffness from tangent stiffness */

Keff = stiff + (2/dt)*damp + (4/dt/dt)*mass;

/* Solve for d_displacement, d_velocity */

dp_i = Solve( Keff, dPeff);
dp = dp + dp_i;
dv = (2/dt)*dp - 2*velocity;

/* Compute displacement, velocity */

new_displ    = displ + dp;
new_velocity = velocity + dv;

/* Compute incremental displacement and internal load using old stiffness */

dFs = stiff*dp_i;

Fs_i = Fs_i + dFs;
if( ii==1 ) {
x = L2Norm(dFs);
if( x == 0 ) { x = 1; }
}

/* Check material yielding and compute new stiffness */

ElmtStateDet( dp_i );
stiff = Stiff();

/* Compute new internal load, damping force, and acceleration */

Fd = damp*new_velocity;
new_accel = mass_inv*( P_ext-Fs-Fd );

/* Calculate the unbalance force, and error percentage */

P_err = Fs_i - Fs;
y     = L2Norm(P_err);
err   = y/x;

/* Assign new effective incremental load */

dPeff   = P_err;
displ_i = new_displ;
Fs_i    = Fs;

print "*** In While Loop" ,ii, ", dFs = ", L2Norm(dFs), ", P_err =" , y , ", err =" ,err, "\n";
ii = ii+1;
if( ii > 50 ) { flag = 1;
err = tol;
}
}

print "*** FINISHED CONVERGENCE AT STEP", stepno, "\n";

/* Calculate the energy balance, Wint(n+1) + T(n+1) = Wext(n+1) */
/*                                                              */
/* [1] column 1 is internal work by stiffness                   */
/* [2] column 2 is internal work by damping                     */
/* [3] column 3 is kinetic energy                               */
/* [4] column 4 is external work                                */

trans_vel     = Trans(velocity);         /* V(i-1) */
trans_vel_new = Trans(new_velocity);     /* V(i)   */
trans_dis_new = Trans(new_displ);        /* D(i)   */

if( stepno == 1 ) then {
Wint_fs = dt/2*(trans_vel_new*Fs);
Wint_fd = dt/2*(trans_vel_new*Fd);
T    = (trans_vel_new*mass*new_velocity)/2;
Wext = dt/2*(trans_vel_new*P_ext);
energy[1][1] = Wint_fs[1][1];
energy[1][2] = Wint_fd[1][1];
energy[1][3] = T[1][1];
energy[1][4] = Wext[1][1];
} else {
Wint_fs = dt/2*(trans_vel*Fs_old + trans_vel_new*Fs);
Wint_fd = dt/2*(trans_vel*Fd_old + trans_vel_new*Fd);
T    = (trans_vel_new*mass*new_velocity)/2;
Wext = dt/2*(trans_vel*P_old + trans_vel_new*P_ext);
energy[stepno][1] = energy[stepno-1][1] + Wint_fs[1][1];
energy[stepno][2] = energy[stepno-1][2] + Wint_fd[1][1];
energy[stepno][3] = T[1][1];
energy[stepno][4] = energy[stepno-1][4] + Wext[1][1];
}

print "*** flag 1\n";

/* Tolerance is satisfied, update histories for this time step */

UpdateResponse();
P_old    = P_ext;
Fs_old   = Fs;
Fd_old   = Fd;
displ    = new_displ;
velocity = new_velocity;
accel    = new_accel;

print "*** flag 2\n";

/* Store node displacement for this time step */

i = 1;
for( node_no=1 ; node_no<=5 ; node_no=node_no+1 ) {
dof = GetDof([node_no]);
deck_displ[stepno][i] = displ[ dof[1][accel_dir] ][1];
i = i+1;
}
i = 1;
for( node_no=7 ; node_no<=9 ; node_no=node_no+1 ) {
dof = GetDof([node_no]);
pier_displ[stepno][i] = displ[ dof[1][accel_dir] ][1];
i = i+1;
}

print "*** flag 3\n";

/* Store element shear force for this time step */

/* In the meantime, comment this section out .....
i = 1;
for( elmt_no=5 ; elmt_no<=9 ; elmt_no=elmt_no+1 ) {
stress = GetStress( [elmt_no], displ );
isolator_shear[stepno][i] = stress[1][accel_dir];
i=i+1;
}
*/

print "*** flag 4\n";

/*
stress = GetStress([bot_elmt2],displ);
bottom_shear[stepno][1] = -stress[2][accel_dir];
stress = GetStress([bot_elmt3],displ);
bottom_shear[stepno][2] = -stress[2][accel_dir];
stress = GetStress([bot_elmt4],displ);
bottom_shear[stepno][3] = -stress[2][accel_dir];
*/

print "*** flag 5\n";

if( flag == 1 ) {
stepno = total_stepno;
}

print "deck_displ : ";
print  deck_displ[stepno][1],deck_displ[stepno][2],deck_displ[stepno][3],
deck_displ[stepno][4],deck_displ[stepno][5],"\n";
print "pier_displ : ";
print  pier_displ[stepno][1], pier_displ[stepno][2],
pier_displ[stepno][3],"\n";
print "isol_shear : ";
print  isolator_shear[stepno][1], isolator_shear[stepno][2],
isolator_shear[stepno][3], isolator_shear[stepno][4],
isolator_shear[stepno][5],"\n";
print "bottom_shear : ";
print  bottom_shear[stepno][1], bottom_shear[stepno][2],
bottom_shear[stepno][3], "\n";
print "energy       : ";
print  energy[stepno][1], energy[stepno][2],
energy[stepno][3], energy[stepno][4],"\n";
}

PrintMatrix(deck_displ,pier_displ);
PrintMatrix(isolator_shear,bottom_shear);
PrintMatrix(energy);

quit;

```

Developed in March 1998 by Mark Austin