OpenSees » SRC » element » forceBeamColumn
HingeRadauBeamIntegration3d.cpp

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.9 $
00022 // $Date: 2004/06/07 23:21:19 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/element/forceBeamColumn/HingeRadauBeamIntegration3d.cpp,v $
00024 
00025 #include <HingeRadauBeamIntegration3d.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 HingeRadauBeamIntegration3d::HingeRadauBeamIntegration3d(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_HingeRadau3d),
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 HingeRadauBeamIntegration3d::HingeRadauBeamIntegration3d():
00049   BeamIntegration(BEAM_INTEGRATION_TAG_HingeRadau3d),
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 HingeRadauBeamIntegration3d::~HingeRadauBeamIntegration3d()
00056 {
00057   // Nothing to do
00058 }
00059 
00060 void
00061 HingeRadauBeamIntegration3d::getSectionLocations(int numSections, double L,
00062                                                  double *xi)
00063 {
00064   xi[0] = 0.0;
00065   xi[1] = 1.0;
00066   for (int i = 2; i < numSections; i++)
00067     xi[i] = 0.0;
00068 }
00069 
00070 void
00071 HingeRadauBeamIntegration3d::getSectionWeights(int numSections, double L,
00072                                                double *wt)
00073 {
00074   double oneOverL = 1.0/L;
00075 
00076   wt[0] = lpI*oneOverL;
00077   wt[1] = lpJ*oneOverL;
00078   for (int i = 2; i < numSections; i++)
00079     wt[i] = 1.0;
00080 }
00081 
00082 int
00083 HingeRadauBeamIntegration3d::addElasticFlexibility(double L, Matrix &fElastic)
00084 {
00085   double oneOverL = 1.0/L;
00086 
00087   double Le = L-lpI-lpJ;
00088   fElastic(0,0) += Le/(E*A);
00089   fElastic(5,5) += Le/(G*J);
00090 
00091   double x[4];
00092   double w[4];
00093   
00094   static const double eight3 = 8.0/3.0;
00095   static const double oneOverRoot3 = 1.0/sqrt(3.0);
00096 
00097   double oneOverEIz = 1.0/(E*Iz);
00098   double oneOverEIy = 1.0/(E*Iy);
00099   
00100   x[0] = eight3*lpI;
00101   w[0] = 3.0*lpI;
00102 
00103   x[1] = L-eight3*lpJ;
00104   w[1] = 3.0*lpJ;
00105 
00106   Le = L-4.0*(lpI+lpJ);
00107 
00108   x[2] = 4.0*lpI + 0.5*Le*(1.0-oneOverRoot3);
00109   w[2] = 0.5*Le;
00110 
00111   x[3] = 4.0*lpI + 0.5*Le*(1.0+oneOverRoot3);
00112   w[3] = w[2];
00113 
00114   double tmpZ = 0.0;
00115   double tmpY = 0.0;
00116   double xL, xL1, wt;
00117   for (int i = 0; i < 4; i++) {
00118     xL  = x[i]*oneOverL;
00119     xL1 = xL-1.0;
00120     wt = w[i]*oneOverEIz;
00121     fElastic(1,1) += xL1*xL1*wt;
00122     fElastic(2,2) += xL*xL*wt;
00123     tmpZ          += xL*xL1*wt;
00124     wt = w[i]*oneOverEIy;
00125     fElastic(3,3) += xL1*xL1*wt;
00126     fElastic(4,4) += xL*xL*wt;
00127     tmpY          += xL*xL1*wt;
00128   }
00129   fElastic(1,2) += tmpZ;
00130   fElastic(2,1) += tmpZ;
00131   fElastic(3,4) += tmpY;
00132   fElastic(4,3) += tmpY;
00133 
00134   return -1;
00135 }
00136 
00137 void
00138 HingeRadauBeamIntegration3d::addElasticDeformations(ElementalLoad *theLoad,
00139                                                     double loadFactor,
00140                                                     double L, double *v0)
00141 {
00142   // Length of elastic interior
00143   double Le = L-lpI-lpJ;
00144   if (Le == 0.0)
00145     return;
00146 
00147   int type;
00148   const Vector &data = theLoad->getData(type, loadFactor);
00149 
00150   double oneOverL = 1.0/L;
00151 
00152   if (type == LOAD_TAG_Beam3dUniformLoad) {
00153     double wy = data(0)*loadFactor;  // Transverse
00154     double wz = data(1)*loadFactor;  // Transverse
00155     double wx = data(2)*loadFactor;  // Axial
00156 
00157     // Accumulate basic deformations due to uniform load on interior
00158     // Midpoint rule for axial
00159     v0[0] += wx*0.5*(L-lpI+lpJ)/(E*A)*Le;
00160 
00161     // Two point Gauss for bending ... will not be exact when
00162     // hinge lengths are not equal, but this is not a big deal!!!
00163     double x1 = lpI + 0.5*Le*(1.0-1.0/sqrt(3.0));
00164     double x2 = lpI + 0.5*Le*(1.0+1.0/sqrt(3.0));
00165 
00166     double Mz1 = 0.5*wy*x1*(x1-L);
00167     double Mz2 = 0.5*wy*x2*(x2-L);
00168 
00169     double My1 = -0.5*wz*x1*(x1-L);
00170     double My2 = -0.5*wz*x2*(x2-L);
00171 
00172     double Le2EIz = Le/(2*E*Iz);
00173     double Le2EIy = Le/(2*E*Iy);
00174 
00175     double b1, b2;
00176     b1 = x1*oneOverL;
00177     b2 = x2*oneOverL;
00178     v0[2] += Le2EIz*(b1*Mz1+b2*Mz2);
00179     v0[4] += Le2EIy*(b1*My1+b2*My2);
00180 
00181     b1 -= 1.0;
00182     b2 -= 1.0;
00183     v0[1] += Le2EIz*(b1*Mz1+b2*Mz2);
00184     v0[3] += Le2EIy*(b1*My1+b2*My2);
00185   }
00186 
00187   else if (type == LOAD_TAG_Beam3dPointLoad) {
00188     double Py = data(0)*loadFactor;
00189     double Pz = data(1)*loadFactor;
00190     double N  = data(2)*loadFactor;
00191     double aOverL = data(3);
00192     double a = aOverL*L;
00193 
00194     double Vy2 = Py*aOverL;
00195     double Vy1 = Py-Vy2;
00196 
00197     double Vz2 = Pz*aOverL;
00198     double Vz1 = Pz-Vz2;
00199 
00200     // Accumulate basic deformations of interior due to point load
00201     double M1, M2, M3;
00202     double b1, b2, b3;
00203 
00204     // Point load is on left hinge
00205     if (a < lpI) {
00206       M1 = (lpI-L)*Vy2;
00207       M2 = -lpJ*Vy2;
00208 
00209       double Le_6EI;
00210       Le_6EI = Le/(6*E*Iz);
00211 
00212       b1 = lpI*oneOverL;
00213       b2 = 1.0-lpJ*oneOverL;
00214       v0[2] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00215 
00216       b1 -= 1.0;
00217       b2 -= 1.0;
00218       v0[1] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00219 
00220       M1 = (L-lpI)*Vz2;
00221       M2 = lpJ*Vz2;
00222 
00223       Le_6EI = Le/(6*E*Iy);
00224 
00225       b1 = lpI*oneOverL;
00226       b2 = 1.0-lpJ*oneOverL;
00227       v0[4] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00228 
00229       b1 -= 1.0;
00230       b2 -= 1.0;
00231       v0[3] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00232 
00233       // Nothing to do for axial
00234       //v0[0] += 0.0;
00235     }
00236     // Point load is on right hinge
00237     else if (a > L-lpJ) {
00238       M1 = -lpI*Vy1;
00239       M2 = (lpJ-L)*Vy1;
00240 
00241       double Le_6EI;
00242       Le_6EI = Le/(6*E*Iz);
00243 
00244       b1 = lpI*oneOverL;
00245       b2 = 1.0-lpJ*oneOverL;
00246       v0[2] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00247 
00248       b1 -= 1.0;
00249       b2 -= 1.0;
00250       v0[1] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00251 
00252       M1 = lpI*Vz1;
00253       M2 = (L-lpJ)*Vz1;
00254 
00255       Le_6EI = Le/(6*E*Iy);
00256 
00257       b1 = lpI*oneOverL;
00258       b2 = 1.0-lpJ*oneOverL;
00259       v0[4] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00260 
00261       b1 -= 1.0;
00262       b2 -= 1.0;
00263       v0[3] += Le_6EI*(M1*(2*b1+b2)+M2*(b1+2*b2));
00264 
00265       v0[0] += N*Le/(E*A);      
00266     }
00267     // Point load is on elastic interior
00268     else {
00269       M1 = -lpI*Vy1;
00270       M2 = -lpJ*Vy2;
00271       M3 = -a*Vy1;
00272 
00273       double L1_6EI;
00274       double L2_6EI;
00275 
00276       L1_6EI = (a-lpI)/(6*E*Iz);
00277       L2_6EI = (Le-a+lpI)/(6*E*Iz);
00278 
00279       b1 = lpI*oneOverL;
00280       b2 = 1.0-lpJ*oneOverL;
00281       b3 = a*oneOverL;
00282       v0[2] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00283       v0[2] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00284 
00285       b1 -= 1.0;
00286       b2 -= 1.0;
00287       b3 -= 1.0;
00288       v0[1] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00289       v0[1] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00290 
00291       M1 = lpI*Vz1;
00292       M2 = lpJ*Vz2;
00293       M3 = a*Vz1;
00294 
00295       L1_6EI = (a-lpI)/(6*E*Iy);
00296       L2_6EI = (Le-a+lpI)/(6*E*Iy);
00297 
00298       b1 = lpI*oneOverL;
00299       b2 = 1.0-lpJ*oneOverL;
00300       b3 = a*oneOverL;
00301       v0[4] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00302       v0[4] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00303 
00304       b1 -= 1.0;
00305       b2 -= 1.0;
00306       b3 -= 1.0;
00307       v0[3] += L1_6EI*(M1*(2*b1+b3)+M3*(b1+2*b3));
00308       v0[3] += L2_6EI*(M2*(2*b2+b3)+M3*(b2+2*b3));
00309 
00310       v0[0] += N*(a-lpI)/(E*A);
00311     }
00312   }
00313 
00314   return;
00315 }
00316 
00317 double
00318 HingeRadauBeamIntegration3d::getTangentDriftI(double L, double LI,
00319                                               double q2, double q3, bool yAxis)
00320 {
00321   double oneOverL = 1.0/L;
00322 
00323   double betaI = 4*lpI*oneOverL;
00324 
00325   double qq2 = (1-betaI)*q2 - betaI*q3;
00326 
00327   betaI = 8.0/3*lpI*oneOverL;
00328 
00329   double qqq2 = (1-betaI)*q2 - betaI*q3;
00330 
00331   if (LI < lpI)
00332     return 0.0;
00333   else {
00334     double EI = yAxis ? E*Iy : E*Iz;
00335     return (3*lpI)*(LI-8.0/3*lpI)*qqq2/EI +
00336       (LI-4*lpI)/3*(LI-4*lpI)*qq2/EI;
00337   }
00338 }
00339 
00340 double
00341 HingeRadauBeamIntegration3d::getTangentDriftJ(double L, double LI,
00342                                               double q2, double q3, bool yAxis)
00343 {
00344   double oneOverL = 1.0/L;
00345 
00346   double betaJ = 4*lpJ*oneOverL;
00347 
00348   double qq3 = (1-betaJ)*q3 - betaJ*q2;
00349 
00350   betaJ = 8.0/3*lpJ*oneOverL;
00351 
00352   double qqq3 = (1-betaJ)*q3 - betaJ*q2;
00353 
00354   if (LI > L-lpJ)
00355     return 0.0;
00356   else {
00357     double EI = yAxis ? E*Iy : E*Iz;
00358     return (3*lpJ)*(L-LI-8.0/3*lpJ)*qqq3/(EI) +
00359       (L-LI-4*lpJ)/3*(L-LI-4*lpJ)*qq3/(EI);
00360   }
00361 }
00362 
00363 BeamIntegration*
00364 HingeRadauBeamIntegration3d::getCopy(void)
00365 {
00366   return new HingeRadauBeamIntegration3d(E, A, Iz, Iy, G, J, lpI, lpJ);
00367 }
00368 
00369 int
00370 HingeRadauBeamIntegration3d::sendSelf(int cTag, Channel &theChannel)
00371 {
00372   static Vector data(8);
00373 
00374   data(0) = E;
00375   data(1) = A;
00376   data(2) = Iz;
00377   data(3) = Iy;
00378   data(4) = G;
00379   data(5) = J;
00380   data(6) = lpI;
00381   data(7) = lpJ;
00382 
00383   int dbTag = this->getDbTag();
00384 
00385   if (theChannel.sendVector(dbTag, cTag, data) < 0) {
00386     opserr << "HingeRadauBeamIntegration3d::sendSelf() - failed to send Vector data\n";
00387     return -1;
00388   }    
00389 
00390   return 0;
00391 }
00392 
00393 int
00394 HingeRadauBeamIntegration3d::recvSelf(int cTag, Channel &theChannel,
00395                                       FEM_ObjectBroker &theBroker)
00396 {
00397   static Vector data(8);
00398 
00399   int dbTag = this->getDbTag();
00400 
00401   if (theChannel.recvVector(dbTag, cTag, data) < 0)  {
00402     opserr << "HingeRadauBeamIntegration3d::recvSelf() - failed to receive Vector data\n";
00403     return -1;
00404   }
00405   
00406   E   = data(0);
00407   A   = data(1);
00408   Iz  = data(2);
00409   Iy  = data(3);
00410   G   = data(4);
00411   J   = data(5);
00412   lpI = data(6);
00413   lpJ = data(7);
00414 
00415   return 0;
00416 }
00417 
00418 int
00419 HingeRadauBeamIntegration3d::setParameter(const char **argv,
00420                                           int argc, Information &info)
00421 {
00422   if (strcmp(argv[0],"E") == 0) {
00423     info.theType = DoubleType;
00424     return 1;
00425   }
00426   else if (strcmp(argv[0],"A") == 0) {
00427     info.theType = DoubleType;
00428     return 2;
00429   }
00430   else if (strcmp(argv[0],"Iz") == 0) {
00431     info.theType = DoubleType;
00432     return 3;
00433   }
00434   else if (strcmp(argv[0],"Iy") == 0) {
00435     info.theType = DoubleType;
00436     return 4;
00437   }
00438   else if (strcmp(argv[0],"G") == 0) {
00439     info.theType = DoubleType;
00440     return 5;
00441   }
00442   else if (strcmp(argv[0],"J") == 0) {
00443     info.theType = DoubleType;
00444     return 6;
00445   }
00446   else if (strcmp(argv[0],"lpI") == 0) {
00447     info.theType = DoubleType;
00448     return 7;
00449   }
00450   else if (strcmp(argv[0],"lpJ") == 0) {
00451     info.theType = DoubleType;
00452     return 8;
00453   }
00454   else 
00455     return -1;
00456 }
00457 
00458 int
00459 HingeRadauBeamIntegration3d::updateParameter(int parameterID,
00460                                              Information &info)
00461 {
00462   switch (parameterID) {
00463   case 1:
00464     E = info.theDouble;
00465     return 0;
00466   case 2:
00467     A = info.theDouble;
00468     return 0;
00469   case 3:
00470     Iz = info.theDouble;
00471     return 0;
00472   case 4:
00473     Iy = info.theDouble;
00474     return 0;
00475   case 5:
00476     G = info.theDouble;
00477     return 0;
00478   case 6:
00479     J = info.theDouble;
00480     return 0;
00481   case 7:
00482     lpI = info.theDouble;
00483     return 0;
00484   case 8:
00485     lpJ = info.theDouble;
00486     return 0;
00487   default:
00488     return -1;
00489   }
00490 }
00491 
00492 int
00493 HingeRadauBeamIntegration3d::activateParameter(int parameterID)
00494 {
00495   // For Terje to do
00496   return 0;
00497 }
00498 
00499 void
00500 HingeRadauBeamIntegration3d::Print(OPS_Stream &s, int flag)
00501 {
00502   s << "HingeRadau3d" << endln;
00503   s << " E = " << E;
00504   s << " A = " << A;
00505   s << " Iz = " << Iz;
00506   s << " Iy = " << Iy;
00507   s << " G = " << G;
00508   s << " J = " << J << endln;
00509   s << " lpI = " << lpI;
00510   s << " lpJ = " << lpJ << endln;
00511 
00512   return;
00513 }
Generated on Mon Oct 23 15:05:06 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation