OpenSees » SRC » element » forceBeamColumn

HingeRadauTwoBeamIntegration3d.cpp

Go to the documentation of this file.
00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020 
00021 // $Revision: 1.8 $
00022 // $Date: 2004/06/07 23:21:19 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/element/forceBeamColumn/HingeRadauTwoBeamIntegration3d.cpp,v $
00024 
00025 #include <HingeRadauTwoBeamIntegration3d.h>
00026 #include <ElementalLoad.h>
00027 
00028 #include <Matrix.h>
00029 #include <Vector.h>
00030 #include <Channel.h>
00031 #include <FEM_ObjectBroker.h>
00032 #include <Information.h>
00033 
00034 HingeRadauTwoBeamIntegration3d::HingeRadauTwoBeamIntegration3d(double e,
00035                                                          double a,
00036                                                          double iz,
00037                                                          double iy,
00038                                                          double g,
00039                                                          double j,
00040                                                          double lpi,
00041                                                          double lpj):
00042   BeamIntegration(BEAM_INTEGRATION_TAG_HingeRadauTwo3d),
00043   E(e), A(a), Iz(iz), Iy(iy), G(g), J(j), lpI(lpi), lpJ(lpj)
00044 {
00045   // Nothing to do
00046 }
00047 
00048 HingeRadauTwoBeamIntegration3d::HingeRadauTwoBeamIntegration3d():
00049   BeamIntegration(BEAM_INTEGRATION_TAG_HingeRadauTwo3d),
00050   E(0.0), A(0.0), Iz(0.0), Iy(0.0), G(0.0), J(0.0), lpI(0.0), lpJ(0.0)
00051 {
00052 
00053 }
00054 
00055 HingeRadauTwoBeamIntegration3d::~HingeRadauTwoBeamIntegration3d()
00056 {
00057   // Nothing to do
00058 }
00059 
00060 void
00061 HingeRadauTwoBeamIntegration3d::getSectionLocations(int numSections, double L,
00062                                                  double *xi)
00063 {
00064   double two3oneOverL = (2.0/3.0)/L;
00065 
00066   xi[0] = 0.0;
00067   xi[1] = lpI*two3oneOverL;
00068   xi[2] = 1.0-lpJ*two3oneOverL;
00069   xi[3] = 1.0;
00070   for (int i = 4; i < numSections; i++)
00071     xi[i] = 0.0;
00072 }
00073 
00074 void
00075 HingeRadauTwoBeamIntegration3d::getSectionWeights(int numSections, double L,
00076                                                double *wt)
00077 {
00078   double oneOverL = 1.0/L;
00079 
00080   wt[0] = 0.25*lpI*oneOverL;
00081   wt[1] = 3.0*wt[0];
00082   wt[3] = 0.25*lpJ*oneOverL;
00083   wt[2] = 3.0*wt[3];
00084   for (int i = 4; i < numSections; i++)
00085     wt[i] = 1.0;
00086 }
00087 
00088 int
00089 HingeRadauTwoBeamIntegration3d::addElasticFlexibility(double L, Matrix &fElastic)
00090 {
00091   double oneOverL = 1.0/L;
00092 
00093   // Length of elastic interior
00094   double Le = L-lpI-lpJ;
00095   double Lover3EI = Le/(3*E*Iz);
00096   double Lover6EI = 0.5*Lover3EI;
00097   
00098   // Elastic flexibility of element interior
00099   static Matrix fe(2,2);
00100   fe(0,0) = fe(1,1) =  Lover3EI;
00101   fe(0,1) = fe(1,0) = -Lover6EI;
00102   
00103   // Equilibrium transformation matrix
00104   static Matrix B(2,2);
00105   double betaI = lpI*oneOverL;
00106   double betaJ = lpJ*oneOverL;
00107   B(0,0) = 1.0 - betaI;
00108   B(1,1) = 1.0 - betaJ;
00109   B(0,1) = -betaI;
00110   B(1,0) = -betaJ;
00111   
00112   // Transform the elastic flexibility of the element
00113   // interior to the basic system
00114   static Matrix ftmp(2,2);
00115   ftmp.addMatrixTripleProduct(0.0, B, fe, 1.0);
00116 
00117   fElastic(1,1) += ftmp(0,0);
00118   fElastic(1,2) += ftmp(0,1);
00119   fElastic(2,1) += ftmp(1,0);
00120   fElastic(2,2) += ftmp(1,1);
00121 
00122   Lover3EI = Le/(3*E*Iy);
00123   Lover6EI = 0.5*Lover3EI;
00124   fe(0,0) = fe(1,1) =  Lover3EI;
00125   fe(0,1) = fe(1,0) = -Lover6EI;
00126   ftmp.addMatrixTripleProduct(0.0, B, fe, 1.0);
00127   fElastic(3,3) += ftmp(0,0);
00128   fElastic(3,4) += ftmp(0,1);
00129   fElastic(4,3) += ftmp(1,0);
00130   fElastic(4,4) += ftmp(1,1);
00131 
00132   fElastic(0,0) += Le/(E*A);
00133   fElastic(5,5) += Le/(G*J);
00134 
00135   return -1;
00136 }
00137 
00138 void
00139 HingeRadauTwoBeamIntegration3d::addElasticDeformations(ElementalLoad *theLoad,
00140                                                        double loadFactor,
00141                                                        double L, double *v0)
00142 {
00143   // Length of elastic interior
00144   double Le = L-lpI-lpJ;
00145   if (Le == 0.0)
00146     return;
00147 
00148   int type;
00149   const Vector &data = theLoad->getData(type, loadFactor);
00150 
00151   double oneOverL = 1.0/L;
00152 
00153   if (type == LOAD_TAG_Beam3dUniformLoad) {
00154     double wy = data(0)*loadFactor;  // Transverse
00155     double wz = data(1)*loadFactor;  // Transverse
00156     double wx = data(2)*loadFactor;  // Axial
00157 
00158     // Accumulate basic deformations due to uniform load on interior
00159     // Midpoint rule for axial
00160     v0[0] += wx*0.5*(L-lpI+lpJ)/(E*A)*Le;
00161 
00162     // Two point Gauss for bending ... will not be exact when
00163     // hinge lengths are not equal, but this is not a big deal!!!
00164     double x1 = lpI + 0.5*Le*(1.0-1.0/sqrt(3.0));
00165     double x2 = lpI + 0.5*Le*(1.0+1.0/sqrt(3.0));
00166 
00167     double Mz1 = 0.5*wy*x1*(x1-L);
00168     double Mz2 = 0.5*wy*x2*(x2-L);
00169 
00170     double My1 = -0.5*wz*x1*(x1-L);
00171     double My2 = -0.5*wz*x2*(x2-L);
00172 
00173     double Le2EIz = Le/(2*E*Iz);
00174     double Le2EIy = Le/(2*E*Iy);
00175 
00176     double b1, b2;
00177     b1 = x1*oneOverL;
00178     b2 = x2*oneOverL;
00179     v0[2] += Le2EIz*(b1*Mz1+b2*Mz2);
00180     v0[4] += Le2EIy*(b1*My1+b2*My2);
00181 
00182     b1 -= 1.0;
00183     b2 -= 1.0;
00184     v0[1] += Le2EIz*(b1*Mz1+b2*Mz2);
00185     v0[3] += Le2EIy*(b1*My1+b2*My2);
00186   }
00187 
00188   else if (type == LOAD_TAG_Beam3dPointLoad) {
00189     double Py = data(0)*loadFactor;
00190     double Pz = data(1)*loadFactor;
00191     double N  = data(2)*loadFactor;
00192     double aOverL = data(3);
00193     double a = aOverL*L;
00194 
00195     double Vy2 = Py*aOverL;
00196     double Vy1 = Py-Vy2;
00197 
00198     double Vz2 = Pz*aOverL;
00199     double Vz1 = Pz-Vz2;
00200 
00201     // Accumulate basic deformations of interior due to point load
00202     double M1, M2, M3;
00203     double b1, b2, b3;
00204 
00205     // Point load is on left hinge
00206     if (a < lpI) {
00207       M1 = (lpI-L)*Vy2;
00208       M2 = -lpJ*Vy2;
00209 
00210       double Le_6EI;
00211       Le_6EI = Le/(6*E*Iz);
00212 
00213       b1 = lpI*oneOverL;
00214       b2 = 1.0-lpJ*oneOverL;
00215       v0[2] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00216 
00217       b1 -= 1.0;
00218       b2 -= 1.0;
00219       v0[1] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00220 
00221       M1 = (L-lpI)*Vz2;
00222       M2 = lpJ*Vz2;
00223 
00224       Le_6EI = Le/(6*E*Iy);
00225 
00226       b1 = lpI*oneOverL;
00227       b2 = 1.0-lpJ*oneOverL;
00228       v0[4] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00229 
00230       b1 -= 1.0;
00231       b2 -= 1.0;
00232       v0[3] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00233 
00234       // Nothing to do for axial
00235       //v0[0] += 0.0;
00236     }
00237     // Point load is on right hinge
00238     else if (a > L-lpJ) {
00239       M1 = -lpI*Vy1;
00240       M2 = (lpJ-L)*Vy1;
00241 
00242       double Le_6EI;
00243       Le_6EI = Le/(6*E*Iz);
00244 
00245       b1 = lpI*oneOverL;
00246       b2 = 1.0-lpJ*oneOverL;
00247       v0[2] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00248 
00249       b1 -= 1.0;
00250       b2 -= 1.0;
00251       v0[1] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00252 
00253       M1 = lpI*Vz1;
00254       M2 = (L-lpJ)*Vz1;
00255 
00256       Le_6EI = Le/(6*E*Iy);
00257 
00258       b1 = lpI*oneOverL;
00259       b2 = 1.0-lpJ*oneOverL;
00260       v0[4] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00261 
00262       b1 -= 1.0;
00263       b2 -= 1.0;
00264       v0[3] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00265 
00266       v0[0] += N*Le/(E*A);      
00267     }
00268     // Point load is on elastic interior
00269     else {
00270       M1 = -lpI*Vy1;
00271       M2 = -lpJ*Vy2;
00272       M3 = -a*Vy1;
00273 
00274       double L1_6EI;
00275       double L2_6EI;
00276 
00277       L1_6EI = (a-lpI)/(6*E*Iz);
00278       L2_6EI = (Le-a+lpI)/(6*E*Iz);
00279 
00280       b1 = lpI*oneOverL;
00281       b2 = 1.0-lpJ*oneOverL;
00282       b3 = a*oneOverL;
00283       v0[2] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00284       v0[2] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00285 
00286       b1 -= 1.0;
00287       b2 -= 1.0;
00288       b3 -= 1.0;
00289       v0[1] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00290       v0[1] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00291 
00292       M1 = lpI*Vz1;
00293       M2 = lpJ*Vz2;
00294       M3 = a*Vz1;
00295 
00296       L1_6EI = (a-lpI)/(6*E*Iy);
00297       L2_6EI = (Le-a+lpI)/(6*E*Iy);
00298 
00299       b1 = lpI*oneOverL;
00300       b2 = 1.0-lpJ*oneOverL;
00301       b3 = a*oneOverL;
00302       v0[4] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00303       v0[4] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00304 
00305       b1 -= 1.0;
00306       b2 -= 1.0;
00307       b3 -= 1.0;
00308       v0[3] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00309       v0[3] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00310 
00311       v0[0] += N*(a-lpI)/(E*A);
00312     }
00313   }
00314 
00315   return;
00316 }
00317 
00318 double
00319 HingeRadauTwoBeamIntegration3d::getTangentDriftI(double L, double LI,
00320                                                  double q2, double q3, bool yAxis)
00321 {
00322   double oneOverL = 1.0/L;
00323 
00324   double betaI = lpI*oneOverL;
00325 
00326   double qq2 = (1-betaI)*q2 - betaI*q3;
00327 
00328   if (LI < lpI)
00329     return 0.0;
00330   else {
00331     double EI = yAxis ? E*Iy : E*Iz;
00332     return (LI-lpI)/3*(LI-lpI)*qq2/(EI);
00333   }
00334 }
00335 
00336 double
00337 HingeRadauTwoBeamIntegration3d::getTangentDriftJ(double L, double LI,
00338                                                  double q2, double q3, bool yAxis)
00339 {
00340   double oneOverL = 1.0/L;
00341 
00342   double betaJ = lpJ*oneOverL;
00343 
00344   double qq3 = (1-betaJ)*q3 - betaJ*q2;
00345 
00346   if (LI > L-lpJ)
00347     return 0.0;
00348   else {
00349     double EI = yAxis ? E*Iy : E*Iz;
00350     return (L-LI-lpJ)/3*(L-LI-lpJ)*qq3/EI;
00351   }
00352 }
00353 
00354 BeamIntegration*
00355 HingeRadauTwoBeamIntegration3d::getCopy(void)
00356 {
00357   return new HingeRadauTwoBeamIntegration3d(E, A, Iz, Iy, G, J, lpI, lpJ);
00358 }
00359 
00360 int
00361 HingeRadauTwoBeamIntegration3d::sendSelf(int cTag, Channel &theChannel)
00362 {
00363   static Vector data(8);
00364 
00365   data(0) = E;
00366   data(1) = A;
00367   data(2) = Iz;
00368   data(3) = Iy;
00369   data(4) = G;
00370   data(5) = J;
00371   data(6) = lpI;
00372   data(7) = lpJ;
00373 
00374   int dbTag = this->getDbTag();
00375 
00376   if (theChannel.sendVector(dbTag, cTag, data) < 0) {
00377     opserr << "HingeRadauTwoBeamIntegration3d::sendSelf() - failed to send Vector data\n";
00378     return -1;
00379   }    
00380 
00381   return 0;
00382 }
00383 
00384 int
00385 HingeRadauTwoBeamIntegration3d::recvSelf(int cTag, Channel &theChannel,
00386                                       FEM_ObjectBroker &theBroker)
00387 {
00388   static Vector data(8);
00389 
00390   int dbTag = this->getDbTag();
00391 
00392   if (theChannel.recvVector(dbTag, cTag, data) < 0)  {
00393     opserr << "HingeRadauTwoBeamIntegration3d::recvSelf() - failed to receive Vector data\n";
00394     return -1;
00395   }
00396   
00397   E   = data(0);
00398   A   = data(1);
00399   Iz  = data(2);
00400   Iy  = data(3);
00401   G   = data(4);
00402   J   = data(5);
00403   lpI = data(6);
00404   lpJ = data(7);
00405 
00406   return 0;
00407 }
00408 
00409 int
00410 HingeRadauTwoBeamIntegration3d::setParameter(const char **argv,
00411                                              int argc, Information &info)
00412 {
00413   if (strcmp(argv[0],"E") == 0) {
00414     info.theType = DoubleType;
00415     return 1;
00416   }
00417   else if (strcmp(argv[0],"A") == 0) {
00418     info.theType = DoubleType;
00419     return 2;
00420   }
00421   else if (strcmp(argv[0],"Iz") == 0) {
00422     info.theType = DoubleType;
00423     return 3;
00424   }
00425   else if (strcmp(argv[0],"Iy") == 0) {
00426     info.theType = DoubleType;
00427     return 4;
00428   }
00429   else if (strcmp(argv[0],"G") == 0) {
00430     info.theType = DoubleType;
00431     return 5;
00432   }
00433   else if (strcmp(argv[0],"J") == 0) {
00434     info.theType = DoubleType;
00435     return 6;
00436   }
00437   else if (strcmp(argv[0],"lpI") == 0) {
00438     info.theType = DoubleType;
00439     return 7;
00440   }
00441   else if (strcmp(argv[0],"lpJ") == 0) {
00442     info.theType = DoubleType;
00443     return 8;
00444   }
00445   else 
00446     return -1;
00447 }
00448 
00449 int
00450 HingeRadauTwoBeamIntegration3d::updateParameter(int parameterID,
00451                                                 Information &info)
00452 {
00453   switch (parameterID) {
00454   case 1:
00455     E = info.theDouble;
00456     return 0;
00457   case 2:
00458     A = info.theDouble;
00459     return 0;
00460   case 3:
00461     Iz = info.theDouble;
00462     return 0;
00463   case 4:
00464     Iy = info.theDouble;
00465     return 0;
00466   case 5:
00467     G = info.theDouble;
00468     return 0;
00469   case 6:
00470     J = info.theDouble;
00471     return 0;
00472   case 7:
00473     lpI = info.theDouble;
00474     return 0;
00475   case 8:
00476     lpJ = info.theDouble;
00477     return 0;
00478   default:
00479     return -1;
00480   }
00481 }
00482 
00483 int
00484 HingeRadauTwoBeamIntegration3d::activateParameter(int parameterID)
00485 {
00486   // For Terje to do
00487   return 0;
00488 }
00489 
00490 void
00491 HingeRadauTwoBeamIntegration3d::Print(OPS_Stream &s, int flag)
00492 {
00493   s << "HingeRadauTwo3d" << endln;
00494   s << " E = " << E;
00495   s << " A = " << A;
00496   s << " Iz = " << Iz;
00497   s << " Iy = " << Iy;
00498   s << " G = " << G;
00499   s << " J = " << J << endln;
00500   s << " lpI = " << lpI;
00501   s << " lpJ = " << lpJ << endln;
00502 
00503   return;
00504 }

Generated on Mon Oct 23 15:05:06 2006 for OpenSees by doxygen 1.5.0