/* * =========================================================================== * Use direct stiffness method to compute joint displacements and distribution * of bending moments and axial forces in moment resitatant frame. * =========================================================================== * Copyright (C) 1995-96 by Mark Austin and David Mazzoni. * * This software is provided "as is" without express or implied warranty. * Permission is granted to use this software on any computer system, * and to redistribute it freely, subject to the following restrictions: * * 1. The authors are not responsible for the consequences of use of * this software, even if they arise from defects in the software. * 2. The origin of this software must not be misrepresented, either * by explicit claim or by omission. * 3. Altered versions must be plainly marked as such, and must not * be misrepresented as being the original software. * 4. This notice is to remain intact. * * Written By: Mark Austin June 1995 * =========================================================================== */ #include #include #include "matrix.h" #include "vector.h" #include "miscellaneous.h" #include "defs.h" /* Data structures for frame geometry and elements */ typedef struct { float fXcoord; /* x coordinate */ float fYcoord; /* y coordinate */ } NODE; typedef struct { int iEnd1; /* node no at end 1 of beam element */ int iEnd2; /* node no at end 2 of beam element */ float fE; /* Young's modulus of elasticity */ float fInertia; /* section moment of inertia */ float fArea; /* cross section area */ } ELMT; /* * ================================================================= * Functions that define steel moment resistant frame * * NODE *mrfBuildNodeArray() : build array of nodal coordinates. * ELMT *mrfBuildElmtArray() : build array of structural elements. * MATRIX *mrfDestinationArray() : define destination array. * ================================================================= */ enum { NoNodes = 6, NoElmts = 5, NoDofPerElmt = 6 }; NODE *mrfBuildNodeArray() { NODE *spNode; spNode = (NODE *) safeCalloc( NoNodes, sizeof(NODE), __FILE__, __LINE__ ); spNode[0].fXcoord = 0.0; spNode[0].fYcoord = 0.0; spNode[1].fXcoord = 0.0; spNode[1].fYcoord = 10.0; spNode[2].fXcoord = 8.0; spNode[2].fYcoord = 13.0; spNode[3].fXcoord = 12.0; spNode[3].fYcoord = 13.0; spNode[4].fXcoord = 20.0; spNode[4].fYcoord = 10.0; spNode[5].fXcoord = 20.0; spNode[5].fYcoord = 5.0; return ( spNode ); } ELMT *mrfBuildElmtArray() { ELMT *spElmt; int ii; spElmt = (ELMT *) safeCalloc( NoElmts, sizeof(ELMT), __FILE__, __LINE__ ); for (ii = 1; ii <= NoElmts; ii = ii + 1) { spElmt[ ii-1 ].iEnd1 = ii; spElmt[ ii-1 ].iEnd2 = ii+1; spElmt[ ii-1 ].fE = 2.00e+11; spElmt[ ii-1 ].fInertia = 6.60e-04; spElmt[ ii-1 ].fArea = 4.26e-02; } return ( spElmt ); } MATRIX *mrfDestinationArray() { MATRIX *spDest; spDest = matAllocIndirect("Destination Array", IntegerArray, NoElmts , 6); spDest->uMatrix.ipp[ 0 ][ 0 ] = 0; spDest->uMatrix.ipp[ 0 ][ 3 ] = 1; spDest->uMatrix.ipp[ 0 ][ 1 ] = 0; spDest->uMatrix.ipp[ 0 ][ 4 ] = 2; spDest->uMatrix.ipp[ 0 ][ 2 ] = 0; spDest->uMatrix.ipp[ 0 ][ 5 ] = 3; spDest->uMatrix.ipp[ 1 ][ 0 ] = 1; spDest->uMatrix.ipp[ 1 ][ 3 ] = 4; spDest->uMatrix.ipp[ 1 ][ 1 ] = 2; spDest->uMatrix.ipp[ 1 ][ 4 ] = 5; spDest->uMatrix.ipp[ 1 ][ 2 ] = 3; spDest->uMatrix.ipp[ 1 ][ 5 ] = 6; spDest->uMatrix.ipp[ 2 ][ 0 ] = 4; spDest->uMatrix.ipp[ 2 ][ 3 ] = 7; spDest->uMatrix.ipp[ 2 ][ 1 ] = 5; spDest->uMatrix.ipp[ 2 ][ 4 ] = 8; spDest->uMatrix.ipp[ 2 ][ 2 ] = 6; spDest->uMatrix.ipp[ 2 ][ 5 ] = 9; spDest->uMatrix.ipp[ 3 ][ 0 ] = 7; spDest->uMatrix.ipp[ 3 ][ 3 ] = 10; spDest->uMatrix.ipp[ 3 ][ 1 ] = 8; spDest->uMatrix.ipp[ 3 ][ 4 ] = 11; spDest->uMatrix.ipp[ 3 ][ 2 ] = 9; spDest->uMatrix.ipp[ 3 ][ 5 ] = 12; spDest->uMatrix.ipp[ 4 ][ 0 ] = 10; spDest->uMatrix.ipp[ 4 ][ 3 ] = 0; spDest->uMatrix.ipp[ 4 ][ 1 ] = 11; spDest->uMatrix.ipp[ 4 ][ 4 ] = 0; spDest->uMatrix.ipp[ 4 ][ 2 ] = 12; spDest->uMatrix.ipp[ 4 ][ 5 ] = 0; return ( spDest ); } /* Declarations for structural analysis/finite element functions */ MATRIX * feMemberStiffness( double , double , double , double ); MATRIX * feMemberTransformation( double, double, double ); MATRIX * feElementStiffnessMatrix( double, double, double, double, double, double ); int feSizeOfGlobalStiffness ( MATRIX * ); VECTOR * feSkylineProfile( MATRIX *, int ); /* * ==================================================================================== * main() : Use direct stiffness method to compute joint displacements and distribution * of bending moments and axial forces in moment resitatant frame. * ==================================================================================== */ int main( void ) { MATRIX *spStiff, *spKe, *spK; MATRIX *spTemp1, *spTemp2; MATRIX *spElmtDispl, *spDispl; MATRIX *spA, *spDest, *spLoad; VECTOR *spProfile = (VECTOR *) NULL; NODE *spNode; ELMT *spElmt; float fE, fArea, fInertia, fShearForce; float fDx, fDy, fSin, fCos, fElmtLength; int iRowNo, iColumnNo; int iMin, iMax, iMinRow, iSize; int ii, ij, ik, iNodeEnd1, iNodeEnd2; /* [a] : Build nodal coordinate, element property, and destination arrays */ spNode = mrfBuildNodeArray(); spElmt = mrfBuildElmtArray(); spDest = mrfDestinationArray(); matPrint( spDest ); /* [b] : Compute (and print) skyline profile for Global Stiffness Matrix */ iSize = feSizeOfGlobalStiffness( spDest ); spProfile = feSkylineProfile( spDest, iSize ); vecPrint( spProfile ); /* [c] : Assemble Global Stiffness Matrix */ spStiff = matAllocSkyline("Global Stiffness", DoubleArray, iSize, iSize, spProfile->uVector.ip ); for( ii = 1; ii <= NoElmts; ii = ii + 1) { iNodeEnd1 = spElmt[ii-1].iEnd1; iNodeEnd2 = spElmt[ii-1].iEnd2; fDx = spNode[ iNodeEnd2-1 ].fXcoord - spNode[ iNodeEnd1-1 ].fXcoord; fDy = spNode[ iNodeEnd2-1 ].fYcoord - spNode[ iNodeEnd1-1 ].fYcoord; fElmtLength = sqrt( fDx*fDx + fDy*fDy ); fSin = fDy/fElmtLength; fCos = fDx/fElmtLength; spKe = feElementStiffnessMatrix( spElmt[ii-1].fE, spElmt[ii-1].fArea, spElmt[ii-1].fInertia, fElmtLength, fSin, fCos ); for( ij = 1; ij <= NoDofPerElmt; ij = ij + 1) { iRowNo = spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]; if (iRowNo != 0) { for( ik = ij; ik <= NoDofPerElmt; ik = ik + 1) { iColumnNo = spDest->uMatrix.ipp[ ii-1 ][ ik-1 ]; if (iColumnNo != 0) { iMin = MIN( iRowNo, iColumnNo); iMax = MAX( iRowNo, iColumnNo); spStiff->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] += spKe->uMatrix.dpp[ ij-1 ][ ik-1 ]; } } } } } matPrint ( spStiff ); /* [d] : Setup External Load Matrix */ spLoad = matAllocIndirect("External Loads", DoubleArray, iSize, 1 ); spLoad->uMatrix.dpp[ 0 ][ 0 ] = 5000.0; spLoad->uMatrix.dpp[ 1 ][ 0 ] = -4000.0; spLoad->uMatrix.dpp[ 4 ][ 0 ] = -8000.0; spLoad->uMatrix.dpp[ 7 ][ 0 ] = -8000.0; spLoad->uMatrix.dpp[ 10 ][ 0 ] = -4000.0; matPrint ( spLoad ); /* [e] : Compute and print structural displacements */ spTemp1 = matCopy (spStiff); spDispl = matSolve ( spStiff , spLoad ); spDispl->cpName = saveString("Displacements", __FILE__, __LINE__ ); matPrint ( spDispl ); /* [f] : Compute and print element forces */ printf("\nElement Forces\n"); printf("==============\n\n"); spElmtDispl = matAllocIndirect("Elmt Displ", DoubleArray, NoDofPerElmt , 1 ); for( ii = 1; ii <= NoElmts; ii = ii + 1) { /* [f.1] : Transfer global displacements to element coord system */ for( ij = 1; ij <= NoDofPerElmt; ij = ij + 1) { spElmtDispl->uMatrix.dpp[ij-1][0] = 0.0; ik = spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]; if(ik != 0) spElmtDispl->uMatrix.dpp[ij-1][0] = spDispl->uMatrix.dpp[ik-1][0]; } /* [f.2] : Rebuild element stiffness matrix in member coordinate system */ iNodeEnd1 = spElmt[ii-1].iEnd1; iNodeEnd2 = spElmt[ii-1].iEnd2; fDx = spNode[ iNodeEnd2-1 ].fXcoord - spNode[ iNodeEnd1-1 ].fXcoord; fDy = spNode[ iNodeEnd2-1 ].fYcoord - spNode[ iNodeEnd1-1 ].fYcoord; fElmtLength = sqrt( fDx*fDx + fDy*fDy ); fSin = fDy/fElmtLength; fCos = fDx/fElmtLength; fE = spElmt[ii-1].fE; fArea = spElmt[ii-1].fArea; fInertia = spElmt[ii-1].fInertia; /* [f.3] : Member stiffness and geometric transformation matrices */ spK = feMemberStiffness( fE, fArea, fInertia, fElmtLength ); spA = feMemberTransformation( fElmtLength, fSin, fCos); /* [f.4] : Compute forces in member configuration */ spTemp1 = matMult (spA, spElmtDispl ); spTemp2 = matMult (spK, spTemp1 ); /* [f.5] : Print summary of member actions */ fShearForce = (spTemp2->uMatrix.dpp[1][0] + spTemp2->uMatrix.dpp[2][0])/fElmtLength; printf("Element %3d : Axial Force = %12.5g\n", ii, spTemp2->uMatrix.dpp[0][0] ); printf(" : Moment End1 = %12.5g\n", spTemp2->uMatrix.dpp[1][0] ); printf(" : Moment End2 = %12.5g\n", spTemp2->uMatrix.dpp[2][0] ); printf(" : Shear Force = %12.5g\n\n", fShearForce ); matFree ( spTemp1 ); matFree ( spTemp2 ); matFree ( spA ); matFree ( spK ); } } /* * ==================================================================== * feSizeOfGlobalStiffness() : Compute size of global stiffness matrix. * * Input : MATRIX * -- destination array. * Output : int iSize -- no columns/rows in global stiffness matrix. * ==================================================================== */ int feSizeOfGlobalStiffness( MATRIX *spDest ) { int ii, ij, iSize = 0; for( ii = 1; ii <= spDest->iNoRows; ii = ii + 1) { for( ij = 1; ij <= spDest->iNoColumns; ij = ij + 1) { iSize = MAX(iSize, spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]); } } return iSize; } /* * ================================================================================= * feSkylineProfile() : Compute size and skyline profile of global stiffness matrix. * * Input : MATRIX * -- destination array. * : int iSize -- size of Global Stiffness Matrix. * Output : VECTOR * -- height of each column in skyline profile. * ================================================================================= */ VECTOR * feSkylineProfile( MATRIX *spDest, int iSize ) { VECTOR * spProfile; int ii, ij, iMinRow, iColumnNo; spProfile = vecAlloc( "Skyline Profile", IntegerArray, iSize); for( ii = 1; ii <= spDest->iNoRows; ii = ii + 1) { iMinRow = MAX(1, spDest->uMatrix.ipp[ ii-1 ][ 0 ]); for( ij = 2; ij <= NoDofPerElmt; ij = ij + 1) { if( spDest->uMatrix.ipp[ ii-1 ][ ij-1 ] != 0) iMinRow = MIN(iMinRow, spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]); } for( ij = 1; ij <= NoDofPerElmt; ij = ij + 1) { iColumnNo = spDest->uMatrix.ipp[ ii-1 ][ ij-1 ]; if((1 + iColumnNo - iMinRow) > spProfile->uVector.ip[ iColumnNo-1 ]) spProfile->uVector.ip[ iColumnNo-1 ] = 1 + iColumnNo - iMinRow; } } return spProfile; } /* * =================================================================================== * feMemberStiffness() : Compute beam-column element stiffness (in member coordinates) * * Input : dE -- Young's modulus of elasticity. * : dArea -- Section area. * : dInertia -- Section moment of inertia. * : dLength -- Element length. * Output : MATRIX * -- (3x3) beam-column stiffness matrix (in member coordinates). * =================================================================================== */ MATRIX *feMemberStiffness( double dE, double dArea, double dInertia, double dLength ) { MATRIX *spK; spK = matAllocIndirect("[K]" , DoubleArray, 3, 3 ); spK->uMatrix.dpp[ 0 ][ 0 ] = dE*dArea/dLength; spK->uMatrix.dpp[ 1 ][ 1 ] = 4*dE*dInertia/dLength; spK->uMatrix.dpp[ 2 ][ 1 ] = 2*dE*dInertia/dLength; spK->uMatrix.dpp[ 1 ][ 2 ] = 2*dE*dInertia/dLength; spK->uMatrix.dpp[ 2 ][ 2 ] = 4*dE*dInertia/dLength; return (spK); } /* * ======================================================================= * feMemberTransformation() : Compute geometric transformation matrix [a]. * * Input : dLength -- Element length. * : dSin -- sin() of element orientation. * : dCos -- cosine() of element orientation. * Output : MATRIX * -- Geometric transformation matrix [a]. * ======================================================================= */ MATRIX *feMemberTransformation( double dLength, double dSin, double dCos) { MATRIX *spA; spA = matAllocIndirect("[A]" , DoubleArray, 3, 6 ); spA->uMatrix.dpp[ 0 ][ 0 ] = -dCos; spA->uMatrix.dpp[ 1 ][ 0 ] = -dSin/dLength; spA->uMatrix.dpp[ 2 ][ 0 ] = -dSin/dLength; spA->uMatrix.dpp[ 0 ][ 1 ] = -dSin; spA->uMatrix.dpp[ 1 ][ 1 ] = dCos/dLength; spA->uMatrix.dpp[ 2 ][ 1 ] = dCos/dLength; spA->uMatrix.dpp[ 0 ][ 2 ] = 0.0; spA->uMatrix.dpp[ 1 ][ 2 ] = 1.0; spA->uMatrix.dpp[ 2 ][ 2 ] = 0.0; spA->uMatrix.dpp[ 0 ][ 3 ] = dCos; spA->uMatrix.dpp[ 1 ][ 3 ] = dSin/dLength; spA->uMatrix.dpp[ 2 ][ 3 ] = dSin/dLength; spA->uMatrix.dpp[ 0 ][ 4 ] = dSin; spA->uMatrix.dpp[ 1 ][ 4 ] = -dCos/dLength; spA->uMatrix.dpp[ 2 ][ 4 ] = -dCos/dLength; spA->uMatrix.dpp[ 0 ][ 5 ] = 0.0; spA->uMatrix.dpp[ 1 ][ 5 ] = 0.0; spA->uMatrix.dpp[ 2 ][ 5 ] = 1.0; return (spA); } /* * ================================================================================= * feElementStiffnessMatrix() : compute (6x6) beam-column stiffness matrix (in joint * coordinates). * * Input : dE -- Young's modulus of elasticity. * : dArea -- Section area. * : dInertia -- Section moment of inertia. * : dLength -- Element length. * : dSin -- sin() of element orientation. * : dCos -- cosine() of element orientation. * Output : MATRIX * -- (6x6) beam-column stiffness matrix (in joint coordinates). * ================================================================================= */ MATRIX *feElementStiffnessMatrix( double dE, double dArea, double dInertia, double dLength, double dSin, double dCos ) { MATRIX *spA; MATRIX *spK; MATRIX *spKe; MATRIX *spTemp1; MATRIX *spTemp2; /* [a] : Compute member stiffness and geometry transformation matrices */ spK = feMemberStiffness( dE, dArea, dInertia, dLength ); spA = feMemberTransformation( dLength, dSin, dCos); /* [b] : Element stiffness [Ke] = [A]^T.[K].[A] in global coordinate system */ spTemp1 = matTranspose (spA); spTemp2 = matMult (spTemp1, spK ); spKe = matMult (spTemp2, spA ); /* [c] : Release working memory before leaving function */ matFree( spK ); matFree( spA ); matFree( spTemp1 ); matFree( spTemp2 ); return ( spKe ); }