Note : If you download this file and get a run-time stack overflow, edit the file code.c and increase the macro
#define NSTACK 1500
The recompiled program should work just fine.
/* * ============================================================ * Modal Analysis of Five Story Steel Moment Frame subject to * El Centro ground motion. * * Written By: Mark Austin July, 1995 * ============================================================ */ /* [a] : Setup problem specific parameters */ NDimension = 2; NDofPerNode = 3; MaxNodesPerElement = 2; StartMesh(); /* [b] : Generate two-dimensional grid of nodes */ node = 0; for( y = 0 ft; y <= 50 ft; y = y + 10 ft ) { for( x = 0 ft; x <= 55 ft; x = x + 20 ft ) { /* [b.1] : adjust column spacing for central bay */ if(x == 40 ft) { x = x - 5 ft; } /* [b.2] : add new node to finite element mesh */ node = node + 1; AddNode(node, [ x, y ] ); } } /* [c] : Attach column elements to nodes */ elmtno = 0; for ( colno = 1; colno <= 4; colno = colno + 1) { for (floorno = 1; floorno <= 5; floorno = floorno + 1) { elmtno = elmtno + 1; end1 = 4*(floorno - 1) + colno; end2 = end1 + 4; AddElmt( elmtno, [ end1 , end2 ], "mrf_element"); } } /* [d] : Attach beam elements to nodes */ for (floorno = 1; floorno <= 5; floorno = floorno + 1) { for ( bayno = 1; bayno <= 3; bayno = bayno + 1) { end1 = 4*floorno + bayno; end2 = end1 + 1; elmtno = elmtno + 1; AddElmt( elmtno, [ end1 , end2 ], "mrf_element"); } } /* [e] : Define section and material properties */ ElementAttr("mrf_element") { type = "FRAME_2D"; section = "mrfsection"; material = "mrfmaterial"; } SectionAttr("mrfsection") { Izz = 1541.9 in^4; Iyy = 486.3 in^4; depth = 12.0 in; width = 12.0 in; area = 47.4 in^2; } MaterialAttr("mrfmaterial") { density = 0.1024E-5 lb/in^3; poisson = 0.25; yield = 36.0 ksi; E = 29000 ksi; } /* [f] : Apply full-fixity to columns at foundation level */ for(nodeno = 1; nodeno <= 4; nodeno = nodeno + 1) { FixNode( nodeno, [ 1, 1, 1 ]); } for(nodeno = 5; nodeno <= 24; nodeno = nodeno + 1) { FixNode( nodeno, [ 0, 1, 1 ]); } LinkNode([ 5, 6, 7, 8 ], [ 1, 0, 0] ); LinkNode([ 9, 10, 11, 12 ], [ 1, 0, 0] ); LinkNode([ 13, 14, 15, 16 ], [ 1, 0, 0] ); LinkNode([ 17, 18, 19, 20 ], [ 1, 0, 0] ); LinkNode([ 21, 22, 23, 24 ], [ 1, 0, 0] ); /* [g] : Compile and Print Finite Element Mesh */ EndMesh(); PrintMesh(); /* [h] : Compute "stiffness" and "mass" matrices */ stiff = Stiff(); Lumped = [1]; mass = Mass( Lumped ); PrintMatrix(mass); dead_load = 10000 kg; for(i = 1 ; i <= 5; i = i + 1) { mass [i][i] = mass[i][i] + 4*dead_load; } PrintMatrix(stiff); PrintMatrix(mass); /* [i] : Compute and print eigenvalues and eigenvectors */ no_eigen = 2; eigen = Eigen(stiff, mass, [ no_eigen ]); eigenvalue = Eigenvalue(eigen); eigenvector = Eigenvector(eigen); for(i = 1; i <= no_eigen; i = i + 1) { print "Mode", i ," : w^2 = ", eigenvalue[i][1]; print " : T = ", 2*PI/sqrt(eigenvalue[i][1]) ,"\n"; } PrintMatrix(eigenvector); /* [j] : Generalized mass and stiffness matrices */ EigenTrans = Trans(eigenvector); Mstar = EigenTrans*mass*eigenvector; Kstar = EigenTrans*stiff*eigenvector; PrintMatrix( Mstar ); PrintMatrix( Kstar ); /* [k] : Setup Rayleigh Damping for Base - Isolated Structure */ rdamping = 0.03; W1 = sqrt ( eigenvalue[1][1]); W2 = sqrt ( eigenvalue[2][1]); A0 = 2*rdamping*W1*W2/(W1 + W2); A1 = 2*rdamping/(W1 + W2); damp = A0*mass + A1*stiff; PrintMatrix( damp ); /* [k] : Define earthquake loadings... */ Elcentro = ColumnUnits( [ 14.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] ); PrintMatrix( Elcentro ); /* [l] : Initialize system displacement, velocity, and load vectors */ displ = ColumnUnits( Matrix([5,1]), [m] ); vel = ColumnUnits( Matrix([5,1]), [m/sec]); eload = ColumnUnits( Matrix([5,1]), [kN]); r = One([5,1]); /* [m] : Initialize modal displacement, velocity, and acc'n vectors */ Mdispl = ColumnUnits( Matrix([ no_eigen,1 ]), [m] ); Mvel = ColumnUnits( Matrix([ no_eigen,1 ]), [m/sec]); Maccel = ColumnUnits( Matrix([ no_eigen,1 ]), [m/sec/sec]); /* * [n] : Allocate Matrix to store five response parameters -- * Col 1 = time (sec); * Col 2 = 1st mode displacement (cm); * Col 3 = 2nd mode displacement (cm); * Col 4 = 1st + 2nd mode displacement (cm); * Col 5 = Total energy (Joules) */ dt = 0.02 sec; nsteps = 600; beta = 0.25; gamma = 0.5; response = ColumnUnits( Matrix([nsteps+1,5]), [sec], [1]); response = ColumnUnits( response, [cm], [2]); response = ColumnUnits( response, [cm], [3]); response = ColumnUnits( response, [cm], [4]); response = ColumnUnits( response, [Jou], [5]); /* [o] : Compute (and compute LU decomposition) effective mass */ MASS = Mstar + Kstar*beta*dt*dt; lu = Decompose(Copy(MASS)); /* [p] : Mode-Displacement Solution for Response of Undamped MDOF System */ MassTemp = -mass*r; for(i = 1; i <= nsteps; i = i + 1) { print "*** Start Step ",i,"\n"; if(i == 2) { SetUnitsOff; } /* [p.1] : Update external load */ if(i <= 500) then { eload = MassTemp*Elcentro[i][1]; } else { eload = MassTemp*(0); } Pstar = EigenTrans*eload; R = Pstar - Kstar*(Mdispl + Mvel*dt + Maccel*(dt*dt/2.0)*(1-2*beta)); /* [p.2] : Compute new acceleration, velocity and displacement */ Maccel_new = Substitution(lu,R); Mvel_new = Mvel + dt*(Maccel*(1.0-gamma) + gamma*Maccel_new); Mdispl_new = Mdispl + dt*Mvel + ((1 - 2*beta)*Maccel + 2*beta*Maccel_new)*dt*dt/2; /* [p.3] : Update and print new response */ Maccel = Maccel_new; Mvel = Mvel_new; Mdispl = Mdispl_new; /* [p.4] : Combine Modes */ displ = eigenvector*Mdispl; vel = eigenvector*Mvel; /* [p.5] : Compute Total System Energy */ a1 = Trans(vel); a2 = Trans(displ); e1 = a1*mass*vel; e2 = a2*stiff*displ; energy = 0.5*(e1 + e2); /* [p.6] : Save components of time-history response */ response[i+1][1] = i*dt; /* Time */ response[i+1][2] = eigenvector[1][1]*Mdispl[1][1]; /* 1st mode displacement */ response[i+1][3] = eigenvector[1][2]*Mdispl[2][1]; /* 2nd mode displacement */ response[i+1][4] = displ[1][1]; /* 1st + 2nd mode displacement */ response[i+1][5] = energy[1][1]; /* System Energy */ } /* [q] : Print response matrix and quit */ SetUnitsOn; PrintMatrix(response); quit;