OpenSees » SRC » element » updatedLagrangianBeamColumn

Inelastic2DYS03.cpp

Go to the documentation of this file.
00001 // Inelastic2DYS03.cpp
00003 
00004 #include "Inelastic2DYS03.h"
00005 
00007 // Construction/Destruction
00009 
00010 Inelastic2DYS03::Inelastic2DYS03(int tag, double a_ten, double a_com,
00011                  double e, double iz_pos, double iz_neg,
00012                  int Nd1, int Nd2,
00013                                  YieldSurface_BC *ysEnd1,  YieldSurface_BC *ysEnd2,
00014                                  int rf_algo, bool islinear, double rho)
00015 
00016                         :InelasticYS2DGNL(tag, Nd1, Nd2,
00017                         ysEnd1, ysEnd2, rf_algo, islinear, rho),
00018                         Atens(a_ten), Acomp(a_com), E(e), IzPos(iz_pos), IzNeg(iz_neg), ndisp(6), ndisp_hist(6)
00019 {
00020         ndisp_hist.Zero();
00021         ndisp.Zero();
00022 }
00023 
00024 Inelastic2DYS03::~Inelastic2DYS03()
00025 {
00026         //does nothing
00027 }
00028 
00029 
00030 void Inelastic2DYS03::getLocalStiff(Matrix &K)
00031 {
00032         double L1,L2,I1,I2,A;
00033         Vector ndisp_inc(6);
00034 
00035     getIncrNaturalDisp(ndisp_inc);
00036         ndisp = ndisp_hist + ndisp_inc;
00037 
00038         // getTrialnaturalDisp(ndisp);
00039         
00040     opserr << ndisp;
00041     // opserr << ndisp_hist;
00042     // opserr << ndisp_inc;
00043     opserr << "\a";
00044 
00045         if(ndisp(2)*ndisp(5) < 0  || fabs(ndisp(2)*ndisp(5)) < 1e-10) { //if single curvature
00046                 L1 = L;
00047                 L2 = 0;
00048                 if(ndisp(2) > 0 || ndisp(5) < 0)
00049                         I1 = I2 = IzNeg;
00050                 else I1 = I2 = IzPos;
00051         } else {                //double curvature
00052                 
00053                 if((fabs(ndisp(2)) + fabs(ndisp(5)) < 1e-10))
00054                         L1 = 0;
00055                 else
00056                         L1 = (fabs(ndisp(2))*L) / (fabs(ndisp(2)) + fabs(ndisp(5)));
00057                         L2 = L - L1;
00058                 if(ndisp(2) > 0) {
00059                         I1 = IzNeg;
00060                         I2 = IzPos;
00061                 } else {
00062                         I1 = IzPos;
00063                         I2 = IzNeg;
00064                 }
00065         }
00066 
00067         opserr << L1 << "  " << L2 << "\n";
00068         
00069         if(ndisp(3) < 0) //element is in compression
00070                 A = Acomp;
00071         else                                     //element is in tension
00072                 A = Atens;
00073 
00074         //some factors used in stiffness matrix
00075         double X1 = 
00076 I2*I2*L1*L1*L1*L1+4*I2*L1*L1*L1*L2*I1+6*I2*L1*L1*L2*L2*I1+4*I2*L1*L2*L2*L2*I1+L2*L2*L2*L2*I1*I1;
00077         double X2 = (I2*I1*(I2*L1*L1+2*L2*I1*L1+L2*L2*I1))/(X1);
00078         double X3 = ((L1*I2+L2*I1)*I2*I1)/(X1);
00079         double X4 = ((I2*L1*L1+2*I2*L1*L2+L2*L2*I1)*I2*I1)/(X1);
00080 
00081         //zeros in stiffness matrix
00082     K(0,1) = K(0,2) = K(0,4) = K(0,5)=0;
00083     K(1,0) = K(1,3) = 0;
00084     K(2,0) = K(2,3) = 0;
00085     K(3,1) = K(3,2) = K(3,4) = K(3,5)=0;
00086     K(4,0) = K(4,3) = 0;
00087     K(5,0) = K(5,3) = 0;
00088 
00089         //axial terms
00090         K(0,0) = K(3,3) = (A*E)/(L); 
00091         K(0,3) = K(3,0) = (-A*E)/(L);
00092 
00093         //shear and moment terms
00094         K(1,1) = K(4,4) = 12*E*X3;
00095         K(1,4) = K(4,1) = -12*E*X3;
00096         K(1,2) = K(2,1) = 6*E*X2;
00097         K(1,5) = K(5,1) = 6*E*X4;
00098         K(2,4) = K(4,2) = -6*E*X2;
00099         K(4,5) = K(5,4) = -6*E*X4;
00100         K(2,2) = (4*E*I2*I1*(I2*L1*L1*L1 + 3*L2*I1*L1*L1 + 3*L2*L2*I1*L1 + 
00101               L2*L2*L2*I1))/X1;
00102         K(5,5) = (4*E*I2*I1*(I2*L1*L1*L1 + 3*I2*L1*L1*L2 + 3*I2*L1*L2*L2 + 
00103               L2*L2*L2*I1))/X1;
00104         K(2,5) = K(5,2) = 2*E*I2*I1*(I2*L1*L1*L1 + 3*I2*L1*L1*L2 + 3*L2*L2*I1*L1 
00105              + L2*L2*L2*I1)/X1;
00106 
00107    opserr << "\nInelastic2DYS03::getLocalStiff(..) = \n" << K;
00108 }//getLocalStiff
00109 
00110 
00111 int Inelastic2DYS03::commitState()
00112 {
00113         // first let the super classes do their stuff
00114     this->InelasticYS2DGNL::commitState();
00115         // now set the commit natural disps
00116         ndisp_hist = ndisp;
00117         return 0;
00118 
00119 }
00120

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