OpenSees » SRC » element » forceBeamColumn
ForceBeamColumn2d.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.24 $
00022 // $Date: 2006/10/02 18:32:02 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/element/forceBeamColumn/ForceBeamColumn2d.cpp,v $
00024 
00025 #include <math.h>
00026 #include <stdlib.h>
00027 #include <string.h>
00028 #include <float.h>
00029 
00030 #include <Information.h>
00031 #include <Parameter.h>
00032 #include <ForceBeamColumn2d.h>
00033 #include <MatrixUtil.h>
00034 #include <Domain.h>
00035 #include <Channel.h>
00036 #include <FEM_ObjectBroker.h>
00037 #include <Renderer.h>
00038 #include <math.h>
00039 
00040 #include <ElementResponse.h>
00041 #include <ElementalLoad.h>
00042 
00043 #define  NDM   2         // dimension of the problem (2d)
00044 #define  NND   3         // number of nodal dof's
00045 #define  NEGD  6         // number of element global dof's
00046 #define  NEBD  3         // number of element dof's in the basic system
00047 
00048 Matrix ForceBeamColumn2d::theMatrix(6,6);
00049 Vector ForceBeamColumn2d::theVector(6);
00050 double ForceBeamColumn2d::workArea[100];
00051 
00052 Vector *ForceBeamColumn2d::vsSubdivide = 0;
00053 Matrix *ForceBeamColumn2d::fsSubdivide = 0;
00054 Vector *ForceBeamColumn2d::SsrSubdivide = 0;
00055 
00056 // constructor:
00057 // invoked by a FEM_ObjectBroker, recvSelf() needs to be invoked on this object.
00058 ForceBeamColumn2d::ForceBeamColumn2d(): 
00059   Element(0,ELE_TAG_ForceBeamColumn2d), connectedExternalNodes(2), 
00060   beamIntegr(0), numSections(0), sections(0), crdTransf(0),
00061   rho(0.0), maxIters(0), tol(0.0),
00062   initialFlag(0),
00063   kv(NEBD,NEBD), Se(NEBD),
00064   kvcommit(NEBD,NEBD), Secommit(NEBD),
00065   fs(0), vs(0), Ssr(0), vscommit(0), numEleLoads(0), Ki(0),
00066   parameterID(0)
00067 {
00068   theNodes[0] = 0;  
00069   theNodes[1] = 0;
00070 
00071   if (vsSubdivide == 0)
00072     vsSubdivide  = new Vector [maxNumSections];
00073   if (fsSubdivide == 0)
00074     fsSubdivide  = new Matrix [maxNumSections];
00075   if (SsrSubdivide == 0)
00076     SsrSubdivide  = new Vector [maxNumSections];
00077   if (!vsSubdivide || !fsSubdivide || !SsrSubdivide) {
00078     opserr << "ForceBeamColumn2d::ForceBeamColumn2d() -- failed to allocate Subdivide arrays";   
00079     exit(-1);
00080   }
00081 }
00082 
00083 // constructor which takes the unique element tag, sections,
00084 // and the node ID's of it's nodal end points. 
00085 // allocates the necessary space needed by each object
00086 ForceBeamColumn2d::ForceBeamColumn2d (int tag, int nodeI, int nodeJ,
00087                                       int numSec, SectionForceDeformation **sec,
00088                                       BeamIntegration &bi,
00089                                       CrdTransf2d &coordTransf, double massDensPerUnitLength,
00090                                       int maxNumIters, double tolerance):
00091   Element(tag,ELE_TAG_ForceBeamColumn2d), connectedExternalNodes(2),
00092   beamIntegr(0), numSections(0), sections(0), crdTransf(0),
00093   rho(massDensPerUnitLength),maxIters(maxNumIters), tol(tolerance), 
00094   initialFlag(0),
00095   kv(NEBD,NEBD), Se(NEBD), 
00096   kvcommit(NEBD,NEBD), Secommit(NEBD),
00097   fs(0), vs(0),Ssr(0), vscommit(0), numEleLoads(0), Ki(0),
00098   parameterID(0)
00099 {
00100   theNodes[0] = 0;
00101   theNodes[1] = 0;
00102 
00103   connectedExternalNodes(0) = nodeI;
00104   connectedExternalNodes(1) = nodeJ;    
00105   
00106   beamIntegr = bi.getCopy();
00107   if (beamIntegr == 0) {
00108     opserr << "Error: ForceBeamColumn2d::ForceBeamColumn2d: could not create copy of beam integration object" << endln;
00109     exit(-1);
00110   }
00111   
00112   // get copy of the transformation object   
00113   crdTransf = coordTransf.getCopy(); 
00114   if (crdTransf == 0) {
00115     opserr << "Error: ForceBeamColumn2d::ForceBeamColumn2d: could not create copy of coordinate transformation object" << endln;
00116     exit(-1);
00117   }
00118 
00119   this->setSectionPointers(numSec, sec);
00120   
00121   if (vsSubdivide == 0)
00122     vsSubdivide  = new Vector [maxNumSections];
00123   if (fsSubdivide == 0)
00124     fsSubdivide  = new Matrix [maxNumSections];
00125   if (SsrSubdivide == 0)
00126     SsrSubdivide  = new Vector [maxNumSections];
00127   if (!vsSubdivide || !fsSubdivide || !SsrSubdivide) {
00128     opserr << "ForceBeamColumn2d::ForceBeamColumn2d() -- failed to allocate Subdivide arrays";   
00129     exit(-1);
00130   }
00131 }
00132 
00133 // ~ForceBeamColumn2d():
00134 //      destructor
00135 //      delete must be invoked on any objects created by the object
00136 ForceBeamColumn2d::~ForceBeamColumn2d()
00137 {
00138   if (sections != 0) {
00139     for (int i=0; i < numSections; i++)
00140       if (sections[i] != 0)
00141         delete sections[i];
00142     delete [] sections;
00143   }
00144   
00145   if (fs != 0) 
00146     delete [] fs;
00147   
00148   if (vs != 0) 
00149     delete [] vs;
00150   
00151   if (Ssr != 0) 
00152     delete [] Ssr;
00153   
00154   if (vscommit != 0) 
00155     delete [] vscommit;
00156   
00157   if (crdTransf != 0)
00158     delete crdTransf;
00159 
00160   if (beamIntegr != 0)
00161     delete beamIntegr;
00162   
00163   if (Ki != 0)
00164     delete Ki;
00165 }
00166 
00167 int
00168 ForceBeamColumn2d::getNumExternalNodes(void) const
00169 {
00170   return 2;
00171 }
00172 
00173 const ID &
00174 ForceBeamColumn2d::getExternalNodes(void) 
00175 {
00176   return connectedExternalNodes;
00177 }
00178 
00179 Node **
00180 ForceBeamColumn2d::getNodePtrs()
00181 {
00182   return theNodes;
00183 }
00184 
00185 int
00186 ForceBeamColumn2d::getNumDOF(void) 
00187 {
00188   return NEGD;
00189 }
00190 
00191 void
00192 ForceBeamColumn2d::setDomain(Domain *theDomain)
00193 {
00194   // check Domain is not null - invoked when object removed from a domain
00195   if (theDomain == 0) {
00196     theNodes[0] = 0;
00197     theNodes[1] = 0;
00198     
00199     opserr << "ForceBeamColumn2d::setDomain:  theDomain = 0 ";
00200     exit(0); 
00201   }
00202 
00203   // get pointers to the nodes
00204   
00205   int Nd1 = connectedExternalNodes(0);  
00206   int Nd2 = connectedExternalNodes(1);
00207   
00208   theNodes[0] = theDomain->getNode(Nd1);
00209   theNodes[1] = theDomain->getNode(Nd2);  
00210   
00211   if (theNodes[0] == 0) {
00212     opserr << "ForceBeamColumn2d::setDomain: Nd1: ";
00213     opserr << Nd1 << "does not exist in model\n";
00214     exit(0);
00215   }
00216   
00217   if (theNodes[1] == 0) {
00218     opserr << "ForceBeamColumn2d::setDomain: Nd2: ";
00219     opserr << Nd2 << "does not exist in model\n";
00220     exit(0);
00221   }
00222   
00223   // call the DomainComponent class method 
00224   this->DomainComponent::setDomain(theDomain);
00225   
00226   // ensure connected nodes have correct number of dof's
00227   int dofNode1 = theNodes[0]->getNumberDOF();
00228   int dofNode2 = theNodes[1]->getNumberDOF();
00229   
00230   if ((dofNode1 != NND) || (dofNode2 != NND)) {
00231     opserr << "ForceBeamColumn2d::setDomain(): Nd2 or Nd1 incorrect dof ";
00232     exit(0);
00233   }
00234    
00235   // initialize the transformation
00236   if (crdTransf->initialize(theNodes[0], theNodes[1])) {
00237     opserr << "ForceBeamColumn2d::setDomain(): Error initializing coordinate transformation";  
00238     exit(0);
00239   }
00240     
00241   // get element length
00242   double L = crdTransf->getInitialLength();
00243   if (L == 0.0) {
00244     opserr << "ForceBeamColumn2d::setDomain(): Zero element length:" << this->getTag();  
00245     exit(0);
00246   }
00247 
00248   if (initialFlag == 0) 
00249     this->initializeSectionHistoryVariables();
00250 }
00251 
00252 int
00253 ForceBeamColumn2d::commitState()
00254 {
00255   int err = 0;
00256   int i = 0;
00257 
00258   // call element commitState to do any base class stuff
00259   if ((err = this->Element::commitState()) != 0) {
00260     opserr << "ForceBeamColumn2d::commitState () - failed in base class";
00261   }    
00262   
00263   do {
00264     vscommit[i] = vs[i];
00265     err = sections[i++]->commitState();
00266     
00267   } while (err == 0 && i < numSections);
00268   
00269   if (err)
00270     return err;
00271   
00272   // commit the transformation between coord. systems
00273   if ((err = crdTransf->commitState()) != 0)
00274     return err;
00275   
00276   // commit the element variables state
00277   kvcommit = kv;
00278   Secommit = Se;
00279   
00280   //   initialFlag = 0;  fmk - commented out, see what happens to Example3.1.tcl if uncommented
00281   //                         - i have not a clue why, ask remo if he ever gets in contact with us again!
00282   
00283   return err;
00284 }
00285 
00286 int ForceBeamColumn2d::revertToLastCommit()
00287 {
00288   int err;
00289   int i = 0;
00290   
00291   do {
00292     vs[i] = vscommit[i];
00293     err = sections[i]->revertToLastCommit();
00294     
00295     sections[i]->setTrialSectionDeformation(vs[i]);      
00296     Ssr[i] = sections[i]->getStressResultant();
00297     fs[i]  = sections[i]->getSectionFlexibility();
00298     
00299     i++;
00300   } while (err == 0 && i < numSections);
00301   
00302   
00303   if (err)
00304     return err;
00305   
00306   // revert the transformation to last commit
00307   if ((err = crdTransf->revertToLastCommit()) != 0)
00308     return err;
00309   
00310   // revert the element state to last commit
00311   Se   = Secommit;
00312   kv   = kvcommit;
00313   
00314   initialFlag = 0;
00315   // this->update();
00316   
00317   return err;
00318 }
00319 
00320 int ForceBeamColumn2d::revertToStart()
00321 {
00322   // revert the sections state to start
00323   int err;
00324   int i = 0;
00325   
00326   do {
00327     fs[i].Zero();
00328     vs[i].Zero();
00329     Ssr[i].Zero();
00330     err = sections[i++]->revertToStart();
00331     
00332   } while (err == 0 && i < numSections);
00333   
00334   if (err)
00335     return err;
00336   
00337   // revert the transformation to start
00338   if ((err = crdTransf->revertToStart()) != 0)
00339     return err;
00340   
00341   // revert the element state to start
00342   Secommit.Zero();
00343   kvcommit.Zero();
00344   
00345   Se.Zero();
00346   kv.Zero();
00347   
00348   initialFlag = 0;
00349   // this->update();
00350   return err;
00351 }
00352 
00353 
00354 const Matrix &
00355 ForceBeamColumn2d::getInitialStiff(void)
00356 {
00357   // check for quick return
00358   if (Ki != 0)
00359     return *Ki;
00360 
00361   /*
00362   else
00363     Ki = new Matrix(this->getTangentStiff());
00364   */
00365 
00366   static Matrix f(NEBD, NEBD);   // element flexibility matrix  
00367   this->getInitialFlexibility(f);
00368 
00369   /*
00370   static Matrix I(NEBD,NEBD);   // an identity matrix for matrix inverse  
00371   I.Zero();
00372   for (int i=0; i<NEBD; i++)
00373     I(i,i) = 1.0;
00374   
00375   // calculate element stiffness matrix
00376   // invert3by3Matrix(f, kv);
00377 
00378   static Matrix kvInit(NEBD, NEBD);
00379   if (f.Solve(I, kvInit) < 0)
00380     opserr << "ForceBeamColumn2d::getInitialStiff() -- could not invert flexibility\n";
00381   */
00382 
00383   static Matrix kvInit(NEBD, NEBD);
00384   f.Invert(kvInit);
00385   Ki = new Matrix(crdTransf->getInitialGlobalStiffMatrix(kvInit));
00386 
00387   return *Ki;
00388 }
00389 
00390 const Matrix &
00391 ForceBeamColumn2d::getTangentStiff(void)
00392 {
00393   crdTransf->update();  // Will remove once we clean up the corotational 2d transformation -- MHS
00394   return crdTransf->getGlobalStiffMatrix(kv, Se);
00395 }
00396     
00397 void
00398 ForceBeamColumn2d::computeReactions(double *p0)
00399 {
00400   int type;
00401   double L = crdTransf->getInitialLength();
00402   
00403   for (int i = 0; i < numEleLoads; i++) {
00404     
00405     const Vector &data = eleLoads[i]->getData(type, 1.0);
00406     
00407     if (type == LOAD_TAG_Beam2dUniformLoad) {
00408       double wa = data(1)*1.0;  // Axial
00409       double wy = data(0)*1.0;  // Transverse
00410       
00411       p0[0] -= wa*L;
00412       double V = 0.5*wy*L;
00413       p0[1] -= V;
00414       p0[2] -= V;
00415     }
00416     else if (type == LOAD_TAG_Beam2dPointLoad) {
00417       double P = data(0)*1.0;
00418       double N = data(1)*1.0;
00419       double aOverL = data(2);
00420       
00421       if (aOverL < 0.0 || aOverL > 1.0)
00422         continue;
00423       
00424       double a = aOverL*L;
00425       
00426       double V1 = P*(1.0-aOverL);
00427       double V2 = P*aOverL;
00428       
00429       p0[0] -= N;
00430       p0[1] -= V1;
00431       p0[2] -= V2;
00432     }
00433   }
00434 }
00435 
00436 void
00437 ForceBeamColumn2d::computeReactionSensitivity(double *dp0dh, int gradNumber)
00438 {
00439   int type;
00440   double L = crdTransf->getInitialLength();
00441   
00442   double dLdh = crdTransf->getdLdh();
00443 
00444   for (int i = 0; i < numEleLoads; i++) {
00445     
00446     const Vector &data = eleLoads[i]->getData(type, 1.0);
00447     const Vector &sens = eleLoads[i]->getSensitivityData(gradNumber);
00448 
00449     if (type == LOAD_TAG_Beam2dUniformLoad) {
00450       double wa = data(1)*1.0;  // Axial
00451       double dwadh = sens(0);
00452       double wy = data(0)*1.0;  // Transverse
00453       double dwydh = sens(1);
00454       
00455       //p0[0] -= wa*L;
00456       dp0dh[0] -= wa*dLdh + dwadh*L;
00457 
00458       //double V = 0.5*wy*L;
00459       //p0[1] -= V;
00460       //p0[2] -= V;
00461       double dVdh = 0.5*(wy*dLdh + dwydh*L);
00462       dp0dh[1] -= dVdh;
00463       dp0dh[2] -= dVdh;
00464     }
00465     else if (type == LOAD_TAG_Beam2dPointLoad) {
00466       double P = data(0)*1.0;
00467       double dPdh = sens(0);
00468       double N = data(1)*1.0;
00469       double dNdh = sens(1);
00470       double aOverL = data(2);
00471       double daLdh = sens(2);
00472 
00473       if (aOverL < 0.0 || aOverL > 1.0)
00474         continue;
00475       
00476       //double a = aOverL*L;
00477       
00478       //double V1 = P*(1.0-aOverL);
00479       //double V2 = P*aOverL;
00480       double dV1dh = P*(0.0-daLdh) + dPdh*(1.0-aOverL);
00481       double dV2dh = P*daLdh + dPdh*aOverL;
00482       
00483       //p0[0] -= N;
00484       //p0[1] -= V1;
00485       //p0[2] -= V2;
00486       dp0dh[0] -= dNdh;
00487       dp0dh[1] -= dV1dh;
00488       dp0dh[2] -= dV2dh;
00489     }
00490   }
00491 }
00492 
00493 const Vector &
00494 ForceBeamColumn2d::getResistingForce(void)
00495 {
00496   // Will remove once we clean up the corotational 2d transformation -- MHS
00497   crdTransf->update();
00498   
00499   double p0[3];
00500   Vector p0Vec(p0, 3);
00501   p0Vec.Zero();
00502 
00503   if (numEleLoads > 0)
00504     this->computeReactions(p0);
00505 
00506   return crdTransf->getGlobalResistingForce(Se, p0Vec);
00507 }
00508 
00509 void
00510 ForceBeamColumn2d::initializeSectionHistoryVariables (void)
00511 {
00512   for (int i = 0; i < numSections; i++) {
00513     int order = sections[i]->getOrder();
00514     
00515     fs[i] = Matrix(order,order);
00516     vs[i] = Vector(order);
00517     Ssr[i] = Vector(order);
00518     
00519     vscommit[i] = Vector(order);
00520   }
00521 }
00522 
00523 /********* NEWTON , SUBDIVIDE AND INITIAL ITERATIONS ********************
00524  */
00525 int
00526 ForceBeamColumn2d::update()
00527 {
00528   // if have completed a recvSelf() - do a revertToLastCommit
00529   // to get Ssr, etc. set correctly
00530   if (initialFlag == 2)
00531     this->revertToLastCommit();
00532 
00533   // update the transformation
00534   crdTransf->update();
00535        
00536   // get basic displacements and increments
00537   const Vector &v = crdTransf->getBasicTrialDisp();    
00538 
00539   static Vector dv(NEBD);
00540   dv = crdTransf->getBasicIncrDeltaDisp();    
00541 
00542   if (initialFlag != 0 && dv.Norm() <= DBL_EPSILON && numEleLoads == 0)
00543     return 0;
00544 
00545   static Vector vin(NEBD);
00546   vin = v;
00547   vin -= dv;
00548 
00549   double L = crdTransf->getInitialLength();
00550   double oneOverL  = 1.0/L;  
00551 
00552   double xi[maxNumSections];
00553   beamIntegr->getSectionLocations(numSections, L, xi);
00554   
00555   double wt[maxNumSections];
00556   beamIntegr->getSectionWeights(numSections, L, wt);
00557 
00558   static Vector vr(NEBD);       // element residual displacements
00559   static Matrix f(NEBD,NEBD);   // element flexibility matrix
00560   
00561   static Matrix I(NEBD,NEBD);   // an identity matrix for matrix inverse
00562   double dW;                    // section strain energy (work) norm 
00563   int i, j;
00564   
00565   I.Zero();
00566   for (i=0; i<NEBD; i++)
00567     I(i,i) = 1.0;
00568 
00569   int numSubdivide = 1;
00570   bool converged = false;
00571   static Vector dSe(NEBD);
00572   static Vector dvToDo(NEBD);
00573   static Vector dvTrial(NEBD);
00574   static Vector SeTrial(NEBD);
00575   static Matrix kvTrial(NEBD, NEBD);
00576 
00577   dvToDo = dv;
00578   dvTrial = dvToDo;
00579 
00580   static double factor = 10;
00581 
00582   maxSubdivisions = 4;
00583 
00584   // fmk - modification to get compatable ele forces and deformations 
00585   //   for a change in deformation dV we try first a newton iteration, if
00586   //   that fails we try an initial flexibility iteration on first iteration 
00587   //   and then regular newton, if that fails we use the initial flexiblity
00588   //   for all iterations.
00589   //
00590   //   if they both fail we subdivide dV & try to get compatable forces
00591   //   and deformations. if they work and we have subdivided we apply
00592   //   the remaining dV.
00593 
00594   while (converged == false && numSubdivide <= maxSubdivisions) {
00595    
00596     // try regular newton (if l==0), or
00597     // initial tangent on first iteration then regular newton (if l==1), or 
00598     // initial tangent iterations (if l==2)
00599 
00600     for (int l=0; l<3; l++) {
00601 
00602       //      if (l == 1) l = 2;
00603       SeTrial = Se;
00604       kvTrial = kv;
00605       for (i=0; i<numSections; i++) {
00606         vsSubdivide[i] = vs[i];
00607         fsSubdivide[i] = fs[i];
00608         SsrSubdivide[i] = Ssr[i];
00609       }
00610 
00611       // calculate nodal force increments and update nodal forces      
00612       // dSe = kv * dv;
00613       dSe.addMatrixVector(0.0, kvTrial, dvTrial, 1.0);
00614       SeTrial += dSe;
00615 
00616       if (initialFlag != 2) {
00617 
00618         int numIters = maxIters;
00619         if (l == 1) 
00620           numIters = 10*maxIters; // allow 10 times more iterations for initial tangent
00621         
00622         for (j=0; j <numIters; j++) {
00623           // initialize f and vr for integration
00624           f.Zero();
00625           vr.Zero();
00626 
00627           if (beamIntegr->addElasticFlexibility(L, f) < 0) {
00628             vr(0) += f(0,0)*SeTrial(0);
00629             vr(1) += f(1,1)*SeTrial(1) + f(1,2)*SeTrial(2);
00630             vr(2) += f(2,1)*SeTrial(1) + f(2,2)*SeTrial(2);
00631           }
00632 
00633           double v0[3];
00634           v0[0] = 0.0; v0[1] = 0.0; v0[2] = 0.0;
00635 
00636           for (int ie = 0; ie < numEleLoads; ie++)
00637             beamIntegr->addElasticDeformations(eleLoads[ie], 1.0, L, v0);
00638 
00639           // Add effects of element loads
00640           vr(0) += v0[0];
00641           vr(1) += v0[1];
00642           vr(2) += v0[2];
00643           
00644           for (i=0; i<numSections; i++) {
00645             
00646             int order      = sections[i]->getOrder();
00647             const ID &code = sections[i]->getType();
00648             
00649             static Vector Ss;
00650             static Vector dSs;
00651             static Vector dvs;
00652             static Matrix fb;
00653             
00654             Ss.setData(workArea, order);
00655             dSs.setData(&workArea[order], order);
00656             dvs.setData(&workArea[2*order], order);
00657             fb.setData(&workArea[3*order], order, NEBD);
00658             
00659             double xL  = xi[i];
00660             double xL1 = xL-1.0;
00661             double wtL = wt[i]*L;
00662 
00663             // calculate total section forces
00664             // Ss = b*Se + bp*currDistrLoad;
00665             // Ss.addMatrixVector(0.0, b[i], Se, 1.0);
00666             int ii;
00667             for (ii = 0; ii < order; ii++) {
00668               switch(code(ii)) {
00669               case SECTION_RESPONSE_P:
00670                 Ss(ii) = SeTrial(0);
00671                 break;
00672               case SECTION_RESPONSE_MZ:
00673                 Ss(ii) =  xL1*SeTrial(1) + xL*SeTrial(2);
00674                 break;
00675               case SECTION_RESPONSE_VY:
00676                 Ss(ii) = oneOverL*(SeTrial(1)+SeTrial(2));
00677                 break;
00678               default:
00679                 Ss(ii) = 0.0;
00680                 break;
00681               }
00682             }
00683             
00684             // Add the effects of element loads, if present
00685             // s = b*q + sp
00686             if (numEleLoads > 0)
00687               this->computeSectionForces(Ss, i);
00688 
00689             // dSs = Ss - Ssr[i];
00690             dSs = Ss;
00691             dSs.addVector(1.0, SsrSubdivide[i], -1.0);
00692             
00693             // compute section deformation increments
00694             if (l == 0) {
00695 
00696               //  regular newton 
00697               //    vs += fs * dSs;     
00698 
00699               dvs.addMatrixVector(0.0, fsSubdivide[i], dSs, 1.0);
00700             } else if (l == 2) {
00701 
00702               //  newton with initial tangent if first iteration
00703               //    vs += fs0 * dSs;     
00704               //  otherwise regular newton 
00705               //    vs += fs * dSs;     
00706 
00707               if (j == 0) {
00708                 const Matrix &fs0 = sections[i]->getInitialFlexibility();
00709 
00710                 dvs.addMatrixVector(0.0, fs0, dSs, 1.0);
00711               } else
00712                 dvs.addMatrixVector(0.0, fsSubdivide[i], dSs, 1.0);
00713 
00714             } else {
00715              
00716               //  newton with initial tangent
00717               //    vs += fs0 * dSs;     
00718             
00719               const Matrix &fs0 = sections[i]->getInitialFlexibility();
00720               dvs.addMatrixVector(0.0, fs0, dSs, 1.0);
00721             }
00722             
00723             // set section deformations
00724             if (initialFlag != 0)
00725               vsSubdivide[i] += dvs;
00726             
00727             sections[i]->setTrialSectionDeformation(vsSubdivide[i]);
00728             
00729             // get section resisting forces
00730             SsrSubdivide[i] = sections[i]->getStressResultant();
00731             
00732             // get section flexibility matrix
00733             fsSubdivide[i] = sections[i]->getSectionFlexibility();
00734             
00735             // calculate section residual deformations
00736             // dvs = fs * (Ss - Ssr);
00737             dSs = Ss;
00738             dSs.addVector(1.0, SsrSubdivide[i], -1.0);  // dSs = Ss - Ssr[i];
00739             
00740             dvs.addMatrixVector(0.0, fsSubdivide[i], dSs, 1.0);
00741             
00742             // integrate element flexibility matrix
00743             // f = f + (b^ fs * b) * wtL;
00744             //f.addMatrixTripleProduct(1.0, b[i], fs[i], wtL);
00745             int jj;
00746             const Matrix &fSec = fsSubdivide[i];
00747             fb.Zero();
00748             double tmp;
00749             for (ii = 0; ii < order; ii++) {
00750               switch(code(ii)) {
00751               case SECTION_RESPONSE_P:
00752                 for (jj = 0; jj < order; jj++)
00753                   fb(jj,0) += fSec(jj,ii)*wtL;
00754                 break;
00755               case SECTION_RESPONSE_MZ:
00756                 for (jj = 0; jj < order; jj++) {
00757                   tmp = fSec(jj,ii)*wtL;
00758                   fb(jj,1) += xL1*tmp;
00759                   fb(jj,2) += xL*tmp;
00760                 }
00761                 break;
00762               case SECTION_RESPONSE_VY:
00763                 for (jj = 0; jj < order; jj++) {
00764                   tmp = oneOverL*fSec(jj,ii)*wtL;
00765                   fb(jj,1) += tmp;
00766                   fb(jj,2) += tmp;
00767                 }
00768                 break;
00769               default:
00770                 break;
00771               }
00772             }
00773             
00774             for (ii = 0; ii < order; ii++) {
00775               switch (code(ii)) {
00776               case SECTION_RESPONSE_P:
00777                 for (jj = 0; jj < NEBD; jj++)
00778                   f(0,jj) += fb(ii,jj);
00779                 break;
00780               case SECTION_RESPONSE_MZ:
00781                 for (jj = 0; jj < NEBD; jj++) {
00782                   tmp = fb(ii,jj);
00783                   f(1,jj) += xL1*tmp;
00784                   f(2,jj) += xL*tmp;
00785                 }
00786                 break;
00787               case SECTION_RESPONSE_VY:
00788                 for (jj = 0; jj < NEBD; jj++) {
00789                   tmp = oneOverL*fb(ii,jj);
00790                   f(1,jj) += tmp;
00791                   f(2,jj) += tmp;
00792                 }
00793                 break;
00794               default:
00795                 break;
00796               }
00797             }
00798             
00799             // integrate residual deformations
00800             // vr += (b^ (vs + dvs)) * wtL;
00801             //vr.addMatrixTransposeVector(1.0, b[i], vs[i] + dvs, wtL);
00802             dvs.addVector(1.0, vsSubdivide[i], 1.0);
00803             double dei;
00804             for (ii = 0; ii < order; ii++) {
00805               dei = dvs(ii)*wtL;
00806               switch(code(ii)) {
00807               case SECTION_RESPONSE_P:
00808                 vr(0) += dei;
00809                 break;
00810               case SECTION_RESPONSE_MZ:
00811                 vr(1) += xL1*dei; vr(2) += xL*dei;
00812                 break;
00813               case SECTION_RESPONSE_VY:
00814                 tmp = oneOverL*dei;
00815                 vr(1) += tmp; vr(2) += tmp; 
00816                 break;
00817               default:
00818                 break;
00819               }
00820             }
00821           }
00822           
00823           // calculate element stiffness matrix
00824           // invert3by3Matrix(f, kv);     
00825           if (f.Solve(I, kvTrial) < 0)
00826             opserr << "ForceBeamColumn2d::update() -- could not invert flexibility\n";
00827                                     
00828           // dv = vin + dvTrial  - vr
00829           dv = vin;
00830           dv += dvTrial;
00831           dv -= vr;
00832           
00833           // dv.addVector(1.0, vr, -1.0);
00834           
00835           // dSe = kv * dv;
00836           dSe.addMatrixVector(0.0, kvTrial, dv, 1.0);
00837 
00838           dW = dv ^ dSe; 
00839           
00840           SeTrial += dSe;
00841           
00842           // check for convergence of this interval
00843           if (fabs(dW) < tol) { 
00844             
00845             // set the target displacement
00846             dvToDo -= dvTrial;
00847             vin += dvTrial;
00848             
00849             // check if we have got to where we wanted
00850             if (dvToDo.Norm() <= DBL_EPSILON) {
00851               converged = true;
00852               
00853             } else {  // we convreged but we have more to do
00854               //opserr << dvToDo << dvTrial << endln;
00855               // reset variables for start of next subdivision
00856               dvTrial = dvToDo;
00857               numSubdivide = 1;  // NOTE setting subdivide to 1 again maybe too much
00858             }
00859             
00860             // set kv, vs and Se values
00861             kv = kvTrial;
00862             Se = SeTrial;
00863             
00864             for (int k=0; k<numSections; k++) {
00865               vs[k] = vsSubdivide[k];
00866               fs[k] = fsSubdivide[k];
00867               Ssr[k] = SsrSubdivide[k];
00868             }
00869 
00870             // break out of j & l loops
00871             j = numIters+1;
00872             l = 4;
00873 
00874           } else {   //  if (fabs(dW) < tol) { 
00875 
00876             // if we have failed to convrege for all of our newton schemes
00877             // - reduce step size by the factor specified
00878 
00879             if (j == (numIters-1) && (l == 2)) {
00880               dvTrial /= factor;
00881               numSubdivide++;
00882             }
00883           }
00884         } // for (j=0; j<numIters; j++)
00885       } // if (initialFlag != 2)
00886     } // for (int l=0; l<2; l++)
00887   } // while (converged == false)
00888 
00889 
00890   // if fail to converge we return an error flag & print an error message
00891 
00892   if (converged == false) {
00893     opserr << "WARNING - ForceBeamColumn2d::update - failed to get compatible ";
00894     opserr << "element forces & deformations for element: ";
00895     opserr << this->getTag() << "(dW: << " << dW << ")\n";
00896     return -1;
00897   }
00898 
00899   initialFlag = 1;
00900 
00901   return 0;
00902 }
00903 
00904 void ForceBeamColumn2d::getForceInterpolatMatrix(double xi, Matrix &b, const ID &code)
00905 {
00906   b.Zero();
00907   
00908   double L = crdTransf->getInitialLength();
00909   for (int i = 0; i < code.Size(); i++) {
00910     switch (code(i)) {
00911     case SECTION_RESPONSE_MZ:           // Moment, Mz, interpolation
00912       b(i,1) = xi - 1.0;
00913       b(i,2) = xi;
00914       break;
00915     case SECTION_RESPONSE_P:            // Axial, P, interpolation
00916       b(i,0) = 1.0;
00917       break;
00918     case SECTION_RESPONSE_VY:           // Shear, Vy, interpolation
00919       b(i,1) = b(i,2) = 1.0/L;
00920       break;
00921     default:
00922       break;
00923     }
00924   }
00925 }
00926 
00927 void ForceBeamColumn2d::getDistrLoadInterpolatMatrix(double xi, Matrix &bp, const ID &code)
00928 {
00929   bp.Zero();
00930   
00931   double L = crdTransf->getInitialLength();
00932   for (int i = 0; i < code.Size(); i++) {
00933     switch (code(i)) {
00934     case SECTION_RESPONSE_MZ:           // Moment, Mz, interpolation
00935       bp(i,1) = xi*(xi-1)*L*L/2;
00936       break;
00937     case SECTION_RESPONSE_P:            // Axial, P, interpolation
00938       bp(i,0) = (1-xi)*L;
00939       break;
00940     case SECTION_RESPONSE_VY:           // Shear, Vy, interpolation
00941       bp(i,1) = (xi-0.5)*L;
00942       break;
00943     default:
00944       break;
00945     }
00946   }
00947 }
00948 
00949 const Matrix &
00950 ForceBeamColumn2d::getMass(void)
00951 { 
00952   theMatrix.Zero();
00953   
00954   double L = crdTransf->getInitialLength();
00955   if (rho != 0.0)
00956     theMatrix(0,0) = theMatrix(1,1) = theMatrix(3,3) = theMatrix(4,4) = 0.5*L*rho;
00957   
00958   return theMatrix;
00959 }
00960 
00961 void 
00962 ForceBeamColumn2d::zeroLoad(void)
00963 {
00964   // This is a semi-hack -- MHS
00965   numEleLoads = 0;
00966 
00967   return;
00968 }
00969 
00970 int
00971 ForceBeamColumn2d::addLoad(ElementalLoad *theLoad, double loadFactor)
00972 {
00973   if (numEleLoads < maxNumEleLoads)
00974     eleLoads[numEleLoads++] = theLoad;
00975   else
00976     opserr << "ForceBeamColumn2d::addLoad -- maxNumEleLoads exceeded\n";
00977 
00978   return 0;
00979 }
00980 
00981 void
00982 ForceBeamColumn2d::computeSectionForces(Vector &sp, int isec)
00983 {
00984   int type;
00985 
00986   double L = crdTransf->getInitialLength();
00987 
00988   double xi[maxNumSections];
00989   beamIntegr->getSectionLocations(numSections, L, xi);
00990   double x = xi[isec]*L;
00991 
00992   int order = sections[isec]->getOrder();
00993   const ID &code = sections[isec]->getType();
00994 
00995   for (int i = 0; i < numEleLoads; i++) {
00996 
00997     const Vector &data = eleLoads[i]->getData(type, 1.0);
00998     
00999     if (type == LOAD_TAG_Beam2dUniformLoad) {
01000       double wa = data(1)*1.0;  // Axial
01001       double wy = data(0)*1.0;  // Transverse
01002       
01003       for (int ii = 0; ii < order; ii++) {
01004         
01005         switch(code(ii)) {
01006         case SECTION_RESPONSE_P:
01007           sp(ii) += wa*(L-x);
01008           break;
01009         case SECTION_RESPONSE_MZ:
01010           sp(ii) += wy*0.5*x*(x-L);
01011           break;
01012         case SECTION_RESPONSE_VY:
01013           sp(ii) += wy*(x-0.5*L);
01014           break;
01015         default:
01016           break;
01017         }
01018       }
01019     }
01020     else if (type == LOAD_TAG_Beam2dPointLoad) {
01021       double P = data(0)*1.0;
01022       double N = data(1)*1.0;
01023       double aOverL = data(2);
01024       
01025       if (aOverL < 0.0 || aOverL > 1.0)
01026         continue;
01027       
01028       double a = aOverL*L;
01029       
01030       double V1 = P*(1.0-aOverL);
01031       double V2 = P*aOverL;
01032       
01033       for (int ii = 0; ii < order; ii++) {
01034         
01035         if (x <= a) {
01036           switch(code(ii)) {
01037           case SECTION_RESPONSE_P:
01038             sp(ii) += N;
01039             break;
01040           case SECTION_RESPONSE_MZ:
01041             sp(ii) -= x*V1;
01042             break;
01043           case SECTION_RESPONSE_VY:
01044             sp(ii) -= V1;
01045             break;
01046           default:
01047             break;
01048           }
01049         }
01050         else {
01051           switch(code(ii)) {
01052           case SECTION_RESPONSE_MZ:
01053             sp(ii) -= (L-x)*V2;
01054             break;
01055           case SECTION_RESPONSE_VY:
01056             sp(ii) += V2;
01057             break;
01058           default:
01059             break;
01060           }
01061         }
01062       }
01063     }
01064     else {
01065       opserr << "ForceBeamColumn2d::addLoad -- load type unknown for element with tag: " <<
01066         this->getTag() << endln;
01067     }
01068   }
01069   
01070   // Don't think we need to do this anymore -- MHS
01071   //this->update(); // quick fix -- MHS
01072 }
01073 
01074 void
01075 ForceBeamColumn2d::computeSectionForceSensitivity(Vector &dspdh, int isec,
01076                                                   int gradNumber)
01077 {
01078   int type;
01079 
01080   double L = crdTransf->getInitialLength();
01081   double dLdh = crdTransf->getdLdh();
01082 
01083   double xi[maxNumSections];
01084   beamIntegr->getSectionLocations(numSections, L, xi);
01085 
01086   double dxidh[maxNumSections];
01087   beamIntegr->getLocationsDeriv(numSections, L, dLdh, dxidh);
01088 
01089   double x = xi[isec]*L;
01090   double dxdh = xi[isec]*dLdh + dxidh[isec]*L;
01091       
01092   int order = sections[isec]->getOrder();
01093   const ID &code = sections[isec]->getType();
01094 
01095   for (int i = 0; i < numEleLoads; i++) {
01096 
01097     const Vector &data = eleLoads[i]->getData(type, 1.0);
01098     
01099     if (type == LOAD_TAG_Beam2dUniformLoad) {
01100       double wa = data(1)*1.0;  // Axial
01101       double wy = data(0)*1.0;  // Transverse
01102 
01103       const Vector &sens = eleLoads[i]->getSensitivityData(gradNumber);
01104       double dwadh = sens(1);
01105       double dwydh = sens(0);
01106       
01107       for (int ii = 0; ii < order; ii++) {
01108         
01109         switch(code(ii)) {
01110         case SECTION_RESPONSE_P:
01111           //sp(ii) += wa*(L-x);
01112           dspdh(ii) += dwadh*(L-x) + wa*(dLdh-x) + wa*(L-dxdh);
01113           break;
01114         case SECTION_RESPONSE_MZ:
01115           //sp(ii) += wy*0.5*x*(x-L);
01116           dspdh(ii) += 0.5 * (dwydh*x*(x-L) + wy*dxdh*(x-L) + wy*x*(dxdh-dLdh));
01117           break;
01118         case SECTION_RESPONSE_VY:
01119           //sp(ii) += wy*(x-0.5*L);
01120           dspdh(ii) += dwydh*(x-0.5*L) + wy*(dxdh-0.5*dLdh);
01121           break;
01122         default:
01123           break;
01124         }
01125       }
01126     }
01127     else if (type == LOAD_TAG_Beam2dPointLoad) {
01128       double P = data(0)*1.0;
01129       double N = data(1)*1.0;
01130       double aOverL = data(2);
01131 
01132       const Vector &sens = eleLoads[i]->getSensitivityData(gradNumber);
01133       double dPdh = sens(0);
01134       double dNdh = sens(1);
01135       double daLdh = sens(2);
01136       
01137       if (aOverL < 0.0 || aOverL > 1.0)
01138         continue;
01139       
01140       double a = aOverL*L;
01141 
01142       double V1 = P*(1.0-aOverL);
01143       double V2 = P*aOverL;
01144       double dV1dh = P*(0.0-daLdh) + dPdh*(1.0-aOverL);
01145       double dV2dh = P*daLdh + dPdh*aOverL;
01146 
01147       for (int ii = 0; ii < order; ii++) {
01148         
01149         if (x <= a) {
01150           switch(code(ii)) {
01151           case SECTION_RESPONSE_P:
01152             //sp(ii) += N;
01153             dspdh(ii) += dNdh;
01154             break;
01155           case SECTION_RESPONSE_MZ:
01156             //sp(ii) -= x*V1;
01157             dspdh(ii) -= dxdh*V1 + x*dV1dh;
01158             break;
01159           case SECTION_RESPONSE_VY:
01160             //sp(ii) -= V1;
01161             dspdh(ii) -= dV1dh;
01162             break;
01163           default:
01164             break;
01165           }
01166         }
01167         else {
01168           switch(code(ii)) {
01169           case SECTION_RESPONSE_MZ:
01170             //sp(ii) -= (L-x)*V2;
01171             dspdh(ii) -= (dLdh-dxdh)*V2 + (L-x)*dV2dh;
01172             break;
01173           case SECTION_RESPONSE_VY:
01174             //sp(ii) += V2;
01175             dspdh(ii) -= dV2dh;
01176             break;
01177           default:
01178             break;
01179           }
01180         }
01181       }
01182     }
01183     else {
01184       opserr << "ForceBeamColumn2d::computeSectionForceSensitivity -- load type unknown for element with tag: " <<
01185         this->getTag() << endln;
01186     }
01187   }
01188 }
01189 
01190 int 
01191 ForceBeamColumn2d::addInertiaLoadToUnbalance(const Vector &accel)
01192 {
01193   // Check for a quick return
01194   if (rho == 0.0)
01195     return 0;
01196 
01197   // get R * accel from the nodes
01198   const Vector &Raccel1 = theNodes[0]->getRV(accel);
01199   const Vector &Raccel2 = theNodes[1]->getRV(accel);    
01200 
01201   double L = crdTransf->getInitialLength();
01202   double m = 0.5*rho*L;
01203 
01204   // Should be done through p0[0]
01205   /*
01206   load(0) -= m*Raccel1(0);
01207   load(1) -= m*Raccel1(1);
01208   load(3) -= m*Raccel2(0);
01209   load(4) -= m*Raccel2(1);
01210   */
01211 
01212   return 0;
01213 }
01214 
01215 const Vector &
01216 ForceBeamColumn2d::getResistingForceIncInertia()
01217 {       
01218   // Compute the current resisting force
01219   theVector = this->getResistingForce();
01220 
01221   // Check for a quick return
01222   if (rho != 0.0) {
01223     const Vector &accel1 = theNodes[0]->getTrialAccel();
01224     const Vector &accel2 = theNodes[1]->getTrialAccel();
01225     
01226     double L = crdTransf->getInitialLength();
01227     double m = 0.5*rho*L;
01228     
01229     theVector(0) += m*accel1(0);
01230     theVector(1) += m*accel1(1);
01231     theVector(3) += m*accel2(0);
01232     theVector(4) += m*accel2(1);
01233     
01234     // add the damping forces if rayleigh damping
01235     if (alphaM != 0.0 || betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
01236       theVector += this->getRayleighDampingForces();
01237 
01238   } else {
01239     // add the damping forces if rayleigh damping
01240     if (betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
01241       theVector += this->getRayleighDampingForces();
01242   }
01243 
01244   return theVector;
01245 }
01246 
01247 int
01248 ForceBeamColumn2d::sendSelf(int commitTag, Channel &theChannel)
01249 {  
01250   // place the integer data into an ID
01251   int dbTag = this->getDbTag();
01252   int i, j , k;
01253   int loc = 0;
01254 
01255   static ID idData(11);  // one bigger than needed so no clash later
01256   idData(0) = this->getTag();
01257   idData(1) = connectedExternalNodes(0);
01258   idData(2) = connectedExternalNodes(1);
01259   idData(3) = numSections;
01260   idData(4) = maxIters;
01261   idData(5) = initialFlag;
01262 
01263   idData(6) = crdTransf->getClassTag();
01264   int crdTransfDbTag  = crdTransf->getDbTag();
01265   if (crdTransfDbTag  == 0) {
01266     crdTransfDbTag = theChannel.getDbTag();
01267     if (crdTransfDbTag  != 0) 
01268       crdTransf->setDbTag(crdTransfDbTag);
01269   }
01270   idData(7) = crdTransfDbTag;
01271 
01272   idData(8) = beamIntegr->getClassTag();
01273   int beamIntegrDbTag  = beamIntegr->getDbTag();
01274   if (beamIntegrDbTag  == 0) {
01275     beamIntegrDbTag = theChannel.getDbTag();
01276     if (beamIntegrDbTag  != 0) 
01277       beamIntegr->setDbTag(beamIntegrDbTag);
01278   }
01279   idData(9) = beamIntegrDbTag;
01280 
01281   if (theChannel.sendID(dbTag, commitTag, idData) < 0) {
01282     opserr << "ForceBeamColumn2d::sendSelf() - failed to send ID data\n";
01283     return -1;
01284   }    
01285 
01286   // send the coordinate transformation
01287   
01288   if (crdTransf->sendSelf(commitTag, theChannel) < 0) {
01289     opserr << "ForceBeamColumn2d::sendSelf() - failed to send crdTrans\n";
01290     return -1;
01291   }      
01292 
01293   if (beamIntegr->sendSelf(commitTag, theChannel) < 0) {
01294     opserr << "ForceBeamColumn2d::sendSelf() - failed to send beamIntegr\n";
01295     return -1;
01296   }      
01297 
01298   
01299   //
01300   // send an ID for the sections containing each sections dbTag and classTag
01301   // if section ha no dbTag get one and assign it
01302   //
01303 
01304   ID idSections(2*numSections);
01305   loc = 0;
01306   for (i = 0; i<numSections; i++) {
01307     int sectClassTag = sections[i]->getClassTag();
01308     int sectDbTag = sections[i]->getDbTag();
01309     if (sectDbTag == 0) {
01310       sectDbTag = theChannel.getDbTag();
01311       sections[i]->setDbTag(sectDbTag);
01312     }
01313 
01314     idSections(loc) = sectClassTag;
01315     idSections(loc+1) = sectDbTag;
01316     loc += 2;
01317   }
01318 
01319   if (theChannel.sendID(dbTag, commitTag, idSections) < 0)  {
01320     opserr << "ForceBeamColumn2d::sendSelf() - failed to send ID data\n";
01321     return -1;
01322   }    
01323 
01324   //
01325   // send the sections
01326   //
01327   
01328   for (j = 0; j<numSections; j++) {
01329     if (sections[j]->sendSelf(commitTag, theChannel) < 0) {
01330       opserr << "ForceBeamColumn2d::sendSelf() - section " << 
01331         j << "failed to send itself\n";
01332       return -1;
01333     }
01334   }
01335   
01336   // into a vector place distrLoadCommit, rho, UeCommit, Secommit and kvcommit
01337   int secDefSize = 0;
01338   for (i = 0; i < numSections; i++) {
01339      int size = sections[i]->getOrder();
01340      secDefSize   += size;
01341   }
01342 
01343   Vector dData(1+1+NEBD+NEBD*NEBD+secDefSize); 
01344   loc = 0;
01345 
01346   // place double variables into Vector
01347   dData(loc++) = rho;
01348   dData(loc++) = tol;
01349   
01350   // put  distrLoadCommit into the Vector
01351   //  for (i=0; i<NL; i++) 
01352   //dData(loc++) = distrLoadcommit(i);
01353 
01354   // place kvcommit into vector
01355   for (i=0; i<NEBD; i++) 
01356     dData(loc++) = Secommit(i);
01357 
01358   // place kvcommit into vector
01359   for (i=0; i<NEBD; i++) 
01360      for (j=0; j<NEBD; j++)
01361         dData(loc++) = kvcommit(i,j);
01362   
01363   // place vscommit into vector
01364   for (k=0; k<numSections; k++)
01365      for (i=0; i<sections[k]->getOrder(); i++)
01366         dData(loc++) = (vscommit[k])(i);
01367   
01368   if (theChannel.sendVector(dbTag, commitTag, dData) < 0) {
01369     opserr << "ForceBeamColumn2d::sendSelf() - failed to send Vector data\n";
01370      return -1;
01371   }    
01372 
01373   return 0;
01374 }    
01375 
01376 int
01377 ForceBeamColumn2d::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
01378 {
01379   //
01380   // receive the integer data containing tag, numSections and coord transformation info
01381   //
01382   int dbTag = this->getDbTag();
01383   int i,j,k;
01384   
01385   static ID idData(11); // one bigger than needed 
01386 
01387   if (theChannel.recvID(dbTag, commitTag, idData) < 0)  {
01388     opserr << "ForceBeamColumn2d::recvSelf() - failed to recv ID data\n";
01389     return -1;
01390   }    
01391 
01392   this->setTag(idData(0));
01393   connectedExternalNodes(0) = idData(1);
01394   connectedExternalNodes(1) = idData(2);
01395   maxIters = idData(4);
01396   initialFlag = idData(5);
01397   
01398   int crdTransfClassTag = idData(6);
01399   int crdTransfDbTag = idData(7);
01400 
01401   int beamIntegrClassTag = idData(8);
01402   int beamIntegrDbTag = idData(9);
01403 
01404   // create a new crdTransf object if one needed
01405   if (crdTransf == 0 || crdTransf->getClassTag() != crdTransfClassTag) {
01406       if (crdTransf != 0)
01407           delete crdTransf;
01408 
01409       crdTransf = theBroker.getNewCrdTransf2d(crdTransfClassTag);
01410 
01411       if (crdTransf == 0) {
01412         opserr << "ForceBeamColumn2d::recvSelf() - failed to obtain a CrdTrans object with classTag" <<
01413           crdTransfClassTag << endln;
01414         exit(-1);
01415       }
01416   }
01417 
01418   crdTransf->setDbTag(crdTransfDbTag);
01419 
01420   // invoke recvSelf on the coordTransf object
01421   if (crdTransf->recvSelf(commitTag, theChannel, theBroker) < 0)  
01422   {
01423      opserr << "ForceBeamColumn2d::sendSelf() - failed to recv crdTranf\n";
01424                              
01425      return -3;
01426   }      
01427 
01428   // create a new beamIntegr object if one needed
01429   if (beamIntegr == 0 || beamIntegr->getClassTag() != beamIntegrClassTag) {
01430       if (beamIntegr != 0)
01431           delete beamIntegr;
01432 
01433       beamIntegr = theBroker.getNewBeamIntegration(beamIntegrClassTag);
01434 
01435       if (beamIntegr == 0) {
01436         opserr << "ForceBeamColumn2d::recvSelf() - failed to obtain the beam integration object with classTag" <<
01437           beamIntegrClassTag << endln;
01438         exit(-1);
01439       }
01440   }
01441 
01442   beamIntegr->setDbTag(beamIntegrDbTag);
01443 
01444   // invoke recvSelf on the beamIntegr object
01445   if (beamIntegr->recvSelf(commitTag, theChannel, theBroker) < 0)  
01446   {
01447      opserr << "ForceBeamColumn2d::sendSelf() - failed to recv beam integration\n";
01448                              
01449      return -3;
01450   }      
01451   
01452   //
01453   // recv an ID for the sections containing each sections dbTag and classTag
01454   //
01455 
01456   ID idSections(2*idData(3));
01457   int loc = 0;
01458 
01459   if (theChannel.recvID(dbTag, commitTag, idSections) < 0)  {
01460     opserr << "ForceBeamColumn2d::recvSelf() - failed to recv ID data\n";
01461     return -1;
01462   }    
01463 
01464   //
01465   // now receive the sections
01466   //
01467   if (numSections != idData(3)) {
01468 
01469     //
01470     // we do not have correct number of sections, must delete the old and create
01471     // new ones before can recvSelf on the sections
01472     //
01473 
01474     // delete the old
01475     if (numSections != 0) {
01476       for (int i=0; i<numSections; i++)
01477         delete sections[i];
01478       delete [] sections;
01479     }
01480 
01481     // create a section and recvSelf on it
01482     numSections = idData(3);
01483 
01484     // Delete the old
01485     if (vscommit != 0)
01486       delete [] vscommit;
01487 
01488     // Allocate the right number
01489     vscommit = new Vector[numSections];
01490     if (vscommit == 0) {
01491       opserr << "ForceBeamColumn2d::recvSelf -- failed to allocate vscommit array\n";
01492                               
01493       return -1;
01494     }
01495 
01496     // Delete the old
01497     if (fs != 0)
01498       delete [] fs;
01499 
01500     // Allocate the right number
01501     fs = new Matrix[numSections];  
01502     if (fs == 0) {
01503       opserr << "ForceBeamColumn2d::recvSelf -- failed to allocate fs array\n";
01504       return -1;
01505     }
01506 
01507     // Delete the old
01508     if (vs != 0)
01509       delete [] vs;
01510 
01511     // Allocate the right number
01512     vs = new Vector[numSections];  
01513     if (vs == 0) {
01514       opserr << "ForceBeamColumn2d::recvSelf -- failed to allocate vs array\n";
01515       return -1;
01516     }
01517 
01518     // Delete the old
01519     if (Ssr != 0)
01520       delete [] Ssr;
01521     
01522     // Allocate the right number
01523     Ssr = new Vector[numSections];  
01524     if (Ssr == 0) {
01525       opserr << "ForceBeamColumn2d::recvSelf -- failed to allocate Ssr array\n";
01526       return -1;
01527     }
01528 
01529     // create a new array to hold pointers
01530     sections = new SectionForceDeformation *[idData(3)];
01531     if (sections == 0) {
01532       opserr << "ForceBeamColumn2d::recvSelf() - " << 
01533         "out of memory creating sections array of size" << idData(3) << endln;
01534       exit(-1);
01535     }    
01536 
01537     loc = 0;
01538     
01539     for (i=0; i<numSections; i++) {
01540       int sectClassTag = idSections(loc);
01541       int sectDbTag = idSections(loc+1);
01542       loc += 2;
01543       sections[i] = theBroker.getNewSection(sectClassTag);
01544       if (sections[i] == 0) {
01545         opserr << "ForceBeamColumn2d::recvSelf() - " << 
01546           "Broker could not create Section of class type " << sectClassTag << endln;
01547         exit(-1);
01548       }
01549       sections[i]->setDbTag(sectDbTag);
01550       if (sections[i]->recvSelf(commitTag, theChannel, theBroker) < 0) {
01551         opserr << "ForceBeamColumn2d::recvSelf() - section " << 
01552           i << "failed to recv itself\n";
01553         return -1;
01554       }     
01555     }
01556 
01557     this->initializeSectionHistoryVariables();
01558 
01559   } else {
01560 
01561     // 
01562     // for each existing section, check it is of correct type
01563     // (if not delete old & create a new one) then recvSelf on it
01564     //
01565     
01566     loc = 0;
01567     for (i=0; i<numSections; i++) {
01568       int sectClassTag = idSections(loc);
01569       int sectDbTag = idSections(loc+1);
01570       loc += 2;
01571 
01572       // check of correct type
01573       if (sections[i]->getClassTag() !=  sectClassTag) {
01574         // delete the old section[i] and create a new one
01575         delete sections[i];
01576         sections[i] = theBroker.getNewSection(sectClassTag);
01577         if (sections[i] == 0) {
01578           opserr << "ForceBeamColumn2d::recvSelf() - Broker could not create Section of class type" <<
01579             sectClassTag << endln;
01580           return -1;
01581         }
01582       }
01583 
01584       // recvvSelf on it
01585       sections[i]->setDbTag(sectDbTag);
01586       if (sections[i]->recvSelf(commitTag, theChannel, theBroker) < 0) {
01587         opserr << "ForceBeamColumn2d::recvSelf() - section " << 
01588           i << "failed to recv itself\n";
01589         return -1;
01590       }     
01591     }
01592   }
01593   
01594   // into a vector place distrLoadCommit, rho, UeCommit, Secommit and kvcommit
01595   int secDefSize = 0;
01596   for (int ii = 0; ii < numSections; ii++) {
01597      int size = sections[ii]->getOrder();
01598      secDefSize   += size;
01599   }
01600   
01601   Vector dData(1+1+NEBD+NEBD*NEBD+secDefSize);   
01602   
01603   if (theChannel.recvVector(dbTag, commitTag, dData) < 0)  {
01604     opserr << "ForceBeamColumn2d::sendSelf() - failed to send Vector data\n";
01605     return -1;
01606   }    
01607   
01608   loc = 0;
01609 
01610   // place double variables into Vector
01611   rho = dData(loc++);
01612   tol = dData(loc++);
01613 
01614   // put  distrLoadCommit into the Vector
01615   //for (i=0; i<NL; i++) 
01616   // distrLoad(i) = dData(loc++);
01617 
01618   // place kvcommit into vector
01619   for (i=0; i<NEBD; i++) 
01620     Secommit(i) = dData(loc++);
01621 
01622   // place kvcommit into vector
01623   for (i=0; i<NEBD; i++) 
01624      for (j=0; j<NEBD; j++)
01625         kvcommit(i,j) = dData(loc++);
01626 
01627   kv   = kvcommit;
01628   Se   = Secommit;
01629 
01630   for (k = 0; k < numSections; k++) {
01631     int order = sections[k]->getOrder();
01632 
01633     // place vscommit into vector
01634     vscommit[k] = Vector(order);
01635     for (i = 0; i < order; i++)
01636       (vscommit[k])(i) = dData(loc++);
01637   }
01638 
01639   initialFlag = 2;  
01640 
01641   return 0;
01642 }
01643 
01644 int
01645 ForceBeamColumn2d::getInitialFlexibility(Matrix &fe)
01646 {
01647   fe.Zero();
01648   
01649   double L = crdTransf->getInitialLength();
01650   double oneOverL  = 1.0/L;  
01651   
01652   // Flexibility from elastic interior
01653   beamIntegr->addElasticFlexibility(L, fe);
01654   
01655   double xi[maxNumSections];
01656   beamIntegr->getSectionLocations(numSections, L, xi);
01657   
01658   double wt[maxNumSections];
01659   beamIntegr->getSectionWeights(numSections, L, wt);
01660   
01661   for (int i = 0; i < numSections; i++) {
01662     
01663     int order      = sections[i]->getOrder();
01664     const ID &code = sections[i]->getType();
01665     
01666     Matrix fb(workArea, order, NEBD);
01667     
01668     double xL  = xi[i];
01669     double xL1 = xL-1.0;
01670     double wtL = wt[i]*L;
01671 
01672     const Matrix &fSec = sections[i]->getInitialFlexibility();
01673     fb.Zero();
01674     double tmp;
01675     int ii, jj;
01676     for (ii = 0; ii < order; ii++) {
01677       switch(code(ii)) {
01678       case SECTION_RESPONSE_P:
01679         for (jj = 0; jj < order; jj++)
01680           fb(jj,0) += fSec(jj,ii)*wtL;
01681         break;
01682       case SECTION_RESPONSE_MZ:
01683         for (jj = 0; jj < order; jj++) {
01684           tmp = fSec(jj,ii)*wtL;
01685           fb(jj,1) += xL1*tmp;
01686           fb(jj,2) += xL*tmp;
01687         }
01688         break;
01689       case SECTION_RESPONSE_VY:
01690         for (jj = 0; jj < order; jj++) {
01691           tmp = oneOverL*fSec(jj,ii)*wtL;
01692           fb(jj,1) += tmp;
01693           fb(jj,2) += tmp;
01694         }
01695         break;
01696       default:
01697         break;
01698       }
01699     }
01700     for (ii = 0; ii < order; ii++) {
01701       switch (code(ii)) {
01702       case SECTION_RESPONSE_P:
01703         for (jj = 0; jj < NEBD; jj++)
01704           fe(0,jj) += fb(ii,jj);
01705         break;
01706       case SECTION_RESPONSE_MZ:
01707         for (jj = 0; jj < NEBD; jj++) {
01708           tmp = fb(ii,jj);
01709           fe(1,jj) += xL1*tmp;
01710           fe(2,jj) += xL*tmp;
01711         }
01712         break;
01713       case SECTION_RESPONSE_VY:
01714         for (jj = 0; jj < NEBD; jj++) {
01715           tmp = oneOverL*fb(ii,jj);
01716           fe(1,jj) += tmp;
01717           fe(2,jj) += tmp;
01718         }
01719         break;
01720       default:
01721         break;
01722       }
01723     }
01724   }
01725   
01726   return 0;
01727 }
01728 
01729 void ForceBeamColumn2d::compSectionDisplacements(Vector sectionCoords[], Vector sectionDispls[]) const
01730 {
01731    // get basic displacements and increments
01732    static Vector ub(NEBD);
01733    ub = crdTransf->getBasicTrialDisp();    
01734 
01735    double L = crdTransf->getInitialLength();
01736   
01737    // get integration point positions and weights
01738    //   const Matrix &xi_pt  = quadRule.getIntegrPointCoords(numSections);
01739    // get integration point positions and weights
01740    static double xi_pts[maxNumSections];
01741    beamIntegr->getSectionLocations(numSections, L, xi_pts);
01742 
01743    // setup Vandermode and CBDI influence matrices
01744    int i;
01745    double xi;
01746  
01747    // get CBDI influence matrix
01748    Matrix ls(numSections, numSections);
01749    getCBDIinfluenceMatrix(numSections, xi_pts, L, ls);
01750 
01751    // get section curvatures
01752    Vector kappa(numSections);  // curvature
01753    static Vector vs;              // section deformations 
01754 
01755    for (i=0; i<numSections; i++)
01756    {
01757        // THIS IS VERY INEFFICIENT ... CAN CHANGE LATER
01758        int sectionKey = 0;
01759        const ID &code = sections[i]->getType();
01760        int ii;
01761        for (ii = 0; ii < code.Size(); ii++)
01762            if (code(ii) == SECTION_RESPONSE_MZ)
01763            {
01764                sectionKey = ii;
01765                break;
01766            }
01767 
01768        if (ii == code.Size()) {
01769          opserr << "FATAL NLBeamColumn2d::compSectionDispls - section does not provide Mz response\n";
01770          exit(-1);
01771        }
01772                         
01773        // get section deformations
01774        vs = sections[i]->getSectionDeformation();
01775        kappa(i) = vs(sectionKey);
01776    }
01777 
01778    Vector w(numSections);
01779    static Vector xl(NDM), uxb(NDM);
01780    static Vector xg(NDM), uxg(NDM); 
01781 
01782    // w = ls * kappa;  
01783    w.addMatrixVector (0.0, ls, kappa, 1.0);
01784    
01785    for (i=0; i<numSections; i++)
01786    {
01787       xi = xi_pts[i];
01788 
01789       xl(0) = xi * L;
01790       xl(1) = 0;
01791 
01792       // get section global coordinates
01793       sectionCoords[i] = crdTransf->getPointGlobalCoordFromLocal(xl);
01794 
01795       // compute section displacements
01796       uxb(0) = xi * ub(0); // consider linear variation for axial displacement. CHANGE LATER!!!!!!!!!!
01797       uxb(1) = w(i);
01798              
01799       // get section displacements in global system 
01800       sectionDispls[i] = crdTransf->getPointGlobalDisplFromBasic(xi, uxb);
01801    }           
01802    return;             
01803 }
01804 
01805 void
01806 ForceBeamColumn2d::Print(OPS_Stream &s, int flag)
01807 {
01808   if (flag == 2) {
01809 
01810     s << "#ForceBeamColumn2D\n";
01811 
01812     const Vector &node1Crd = theNodes[0]->getCrds();
01813     const Vector &node2Crd = theNodes[1]->getCrds();    
01814     const Vector &node1Disp = theNodes[0]->getDisp();
01815     const Vector &node2Disp = theNodes[1]->getDisp();    
01816     
01817     s << "#NODE " << node1Crd(0) << " " << node1Crd(1) 
01818       << " " << node1Disp(0) << " " << node1Disp(1) << " " << node1Disp(2) << endln;
01819     
01820     s << "#NODE " << node2Crd(0) << " " << node2Crd(1) 
01821       << " " << node2Disp(0) << " " << node2Disp(1) << " " << node2Disp(2) << endln;
01822 
01823     double P  = Secommit(0);
01824     double M1 = Secommit(1);
01825     double M2 = Secommit(2);
01826     double L = crdTransf->getInitialLength();
01827     double V = (M1+M2)/L;
01828     
01829     double p0[3]; p0[0] = 0.0; p0[1] = 0.0; p0[2] = 0.0;
01830     if (numEleLoads > 0)
01831       this->computeReactions(p0);
01832 
01833     s << "#END_FORCES " << -P+p0[0] << " " << V+p0[1] << " " << M1 << endln;
01834     s << "#END_FORCES " << P << " " << -V+p0[2] << " " << M2 << endln;
01835 
01836     // plastic hinge rotation
01837     static Vector vp(3);
01838     static Matrix fe(3,3);
01839     this->getInitialFlexibility(fe);
01840     vp = crdTransf->getBasicTrialDisp();
01841     vp.addMatrixVector(1.0, fe, Se, -1.0);
01842     s << "#PLASTIC_HINGE_ROTATION " << vp[1] << " " << vp[2] << " " << 0.1*L << " " << 0.1*L << endln;
01843 
01844     // allocate array of vectors to store section coordinates and displacements
01845     static int maxNumSections = 0;
01846     static Vector *coords = 0;
01847     static Vector *displs = 0;
01848     if (maxNumSections < numSections) {
01849       if (coords != 0) 
01850         delete [] coords;
01851       if (displs != 0)
01852         delete [] displs;
01853       
01854       coords = new Vector [numSections];
01855       displs = new Vector [numSections];
01856       
01857       if (!coords) {
01858         opserr << "NLBeamColumn3d::Print() -- failed to allocate coords array";   
01859         exit(-1);
01860       }
01861       
01862       int i;
01863       for (i = 0; i < numSections; i++)
01864         coords[i] = Vector(NDM);
01865       
01866       if (!displs) {
01867         opserr << "NLBeamColumn3d::Print() -- failed to allocate coords array";   
01868         exit(-1);
01869       }
01870       
01871       for (i = 0; i < numSections; i++)
01872         displs[i] = Vector(NDM);
01873       
01874       
01875       maxNumSections = numSections;
01876     }
01877     
01878     // compute section location & displacements
01879     this->compSectionDisplacements(coords, displs);
01880     
01881     // spit out the section location & invoke print on the scetion
01882     for (int i=0; i<numSections; i++) {
01883       s << "#SECTION " << (coords[i])(0) << " " << (coords[i])(1);       
01884       s << " " << (displs[i])(0) << " " << (displs[i])(1) << endln;
01885       sections[i]->Print(s, flag); 
01886     }
01887     
01888   } else {
01889 
01890     s << "\nElement: " << this->getTag() << " Type: ForceBeamColumn2d ";
01891     s << "\tConnected Nodes: " << connectedExternalNodes ;
01892     s << "\tNumber of Sections: " << numSections;
01893     s << "\tMass density: " << rho << endln;
01894     beamIntegr->Print(s, flag);
01895     double P  = Secommit(0);
01896     double M1 = Secommit(1);
01897     double M2 = Secommit(2);
01898     double L = crdTransf->getInitialLength();
01899     double V = (M1+M2)/L;
01900     theVector(1) = V;
01901     theVector(4) = -V;
01902     double p0[3]; p0[0] = 0.0; p0[1] = 0.0; p0[2] = 0.0;
01903     if (numEleLoads > 0)
01904       this->computeReactions(p0);
01905 
01906     s << "\tEnd 1 Forces (P V M): " << -P+p0[0] << " " << V+p0[1] << " " << M1 << endln;
01907     s << "\tEnd 2 Forces (P V M): " << P << " " << -V+p0[2] << " " << M2 << endln;
01908     
01909     if (flag == 1) { 
01910       for (int i = 0; i < numSections; i++)
01911         s << "\numSections "<<i<<" :" << *sections[i];
01912     }
01913   }
01914 }
01915 
01916 OPS_Stream &operator<<(OPS_Stream &s, ForceBeamColumn2d &E)
01917 {
01918   E.Print(s);
01919   return s;
01920 }
01921 
01922 
01923 void
01924 ForceBeamColumn2d::setSectionPointers(int numSec, SectionForceDeformation **secPtrs)
01925 {
01926   if (numSec > maxNumSections) {
01927     opserr << "Error: ForceBeamColumn2d::setSectionPointers -- max number of sections exceeded";
01928   }
01929   
01930   numSections = numSec;
01931   
01932   if (secPtrs == 0) {
01933     opserr << "Error: ForceBeamColumn2d::setSectionPointers -- invalid section pointer";
01934   }       
01935   
01936   sections = new SectionForceDeformation *[numSections];
01937   if (sections == 0) {
01938     opserr << "Error: ForceBeamColumn2d::setSectionPointers -- could not allocate section pointers";
01939   }  
01940   
01941   for (int i = 0; i < numSections; i++) {
01942     
01943     if (secPtrs[i] == 0) {
01944       opserr << "Error: ForceBeamColumn2d::setSectionPointers -- null section pointer " << i << endln;
01945     }
01946     
01947     sections[i] = secPtrs[i]->getCopy();
01948 
01949     if (sections[i] == 0) {
01950       opserr << "Error: ForceBeamColumn2d::setSectionPointers -- could not create copy of section " << i << endln;
01951     }
01952   }
01953   
01954   // allocate section flexibility matrices and section deformation vectors
01955   fs  = new Matrix [numSections];
01956   if (fs == 0) {
01957     opserr << "ForceBeamColumn2d::setSectionPointers -- failed to allocate fs array";
01958   }
01959   
01960   vs = new Vector [numSections];
01961   if (vs == 0) {
01962     opserr << "ForceBeamColumn2d::setSectionPointers -- failed to allocate vs array";
01963   }
01964   
01965   Ssr  = new Vector [numSections];
01966   if (Ssr == 0) {
01967     opserr << "ForceBeamColumn2d::setSectionPointers -- failed to allocate Ssr array";
01968   }
01969   
01970   vscommit = new Vector [numSections];
01971   if (vscommit == 0) {
01972     opserr << "ForceBeamColumn2d::setSectionPointers -- failed to allocate vscommit array";   
01973   }
01974   
01975 }
01976 
01977 int
01978 ForceBeamColumn2d::displaySelf(Renderer &theViewer, int displayMode, float fact)
01979 {
01980   // first determine the end points of the beam based on
01981   // the display factor (a measure of the distorted image)
01982   const Vector &end1Crd = theNodes[0]->getCrds();
01983   const Vector &end2Crd = theNodes[1]->getCrds();       
01984 
01985   static Vector v1(3);
01986   static Vector v2(3);
01987 
01988   if (displayMode >= 0) {
01989     const Vector &end1Disp = theNodes[0]->getDisp();
01990     const Vector &end2Disp = theNodes[1]->getDisp();
01991     
01992     for (int i = 0; i < 2; i++) {
01993       v1(i) = end1Crd(i) + end1Disp(i)*fact;
01994       v2(i) = end2Crd(i) + end2Disp(i)*fact;    
01995     }
01996   } else {
01997     int mode = displayMode  *  -1;
01998     const Matrix &eigen1 = theNodes[0]->getEigenvectors();
01999     const Matrix &eigen2 = theNodes[1]->getEigenvectors();
02000     if (eigen1.noCols() >= mode) {
02001       for (int i = 0; i < 2; i++) {
02002         v1(i) = end1Crd(i) + eigen1(i,mode-1)*fact;
02003         v2(i) = end2Crd(i) + eigen2(i,mode-1)*fact;    
02004       }    
02005     } else {
02006       for (int i = 0; i < 2; i++) {
02007         v1(i) = end1Crd(i);
02008         v2(i) = end2Crd(i);
02009       }    
02010     }
02011   }
02012   
02013   return theViewer.drawLine (v1, v2, 1.0, 1.0);
02014 }
02015 
02016 Response*
02017 ForceBeamColumn2d::setResponse(const char **argv, int argc, Information &eleInformation, OPS_Stream &output)
02018 {
02019   Response *theResponse = 0;
02020 
02021   output.tag("ElementOutput");
02022   output.attr("eleType","ForceBeamColumn2d");
02023   output.attr("eleTag",this->getTag());
02024   output.attr("node1",connectedExternalNodes[0]);
02025   output.attr("node2",connectedExternalNodes[1]);
02026 
02027   // global force - 
02028   if (strcmp(argv[0],"forces") == 0 || strcmp(argv[0],"force") == 0
02029       || strcmp(argv[0],"globalForce") == 0 || strcmp(argv[0],"globalForces") == 0) {
02030 
02031     output.tag("ResponseType","Px_1");
02032     output.tag("ResponseType","Py_1");
02033     output.tag("ResponseType","Mz_1");
02034     output.tag("ResponseType","Px_2");
02035     output.tag("ResponseType","Py_2");
02036     output.tag("ResponseType","Mz_2");
02037 
02038     theResponse =  new ElementResponse(this, 1, theVector);
02039   
02040   
02041   // local force -
02042   } else if (strcmp(argv[0],"localForce") == 0 || strcmp(argv[0],"localForces") == 0) {
02043 
02044     output.tag("ResponseType","N_1");
02045     output.tag("ResponseType","V_1");
02046     output.tag("ResponseType","M_1");
02047     output.tag("ResponseType","N_2");
02048     output.tag("ResponseType","V_2");
02049     output.tag("ResponseType","M_2");
02050 
02051     theResponse =  new ElementResponse(this, 2, theVector);
02052   
02053 
02054   // basic force -
02055   } else if (strcmp(argv[0],"basicForce") == 0 || strcmp(argv[0],"basicForces") == 0) {
02056 
02057     output.tag("ResponseType","N");
02058     output.tag("ResponseType","M_1");
02059     output.tag("ResponseType","M_2");
02060 
02061     theResponse =  new ElementResponse(this, 7, Vector(3));
02062 
02063   // chord rotation -
02064   } else if (strcmp(argv[0],"chordRotation") == 0 || strcmp(argv[0],"chordDeformation") == 0 
02065              || strcmp(argv[0],"basicDeformation") == 0) {
02066 
02067     output.tag("ResponseType","eps");
02068     output.tag("ResponseType","theta_1");
02069     output.tag("ResponseType","theta_2");
02070 
02071     theResponse =  new ElementResponse(this, 3, Vector(3));
02072   
02073   // plastic rotation -
02074   } else if (strcmp(argv[0],"plasticRotation") == 0 || strcmp(argv[0],"plasticDeformation") == 0) {
02075 
02076     output.tag("ResponseType","epsP");
02077     output.tag("ResponseType","thetaP_1");
02078     output.tag("ResponseType","thetaP_2");
02079 
02080     theResponse =  new ElementResponse(this, 4, Vector(3));
02081 
02082   // point of inflection
02083   } else if (strcmp(argv[0],"inflectionPoint") == 0) {
02084     
02085     output.tag("ResponseType","inflectionPoint");
02086 
02087     theResponse =  new ElementResponse(this, 5, 0.0);
02088   
02089   // tangent drift
02090   } else if (strcmp(argv[0],"tangentDrift") == 0) {
02091     theResponse =  new ElementResponse(this, 6, Vector(2));
02092 
02093   // basic forces
02094   } else if (strcmp(argv[0],"basicForce") == 0)
02095     theResponse = new ElementResponse(this, 7, Se);
02096 
02097   /*
02098   // Curvature sensitivity
02099   else if (strcmp(argv[0],"dcurvdh") == 0)
02100     return new ElementResponse(this, 7, Vector(numSections));
02101 
02102   // basic deformation sensitivity
02103   else if (strcmp(argv[0],"dvdh") == 0)
02104     return new ElementResponse(this, 8, Vector(3));
02105   */
02106 
02107   // plastic deformation sensitivity
02108   else if (strcmp(argv[0],"dvpdh") == 0)
02109     return new ElementResponse(this, 9, Vector(3));
02110 
02111   else if (strcmp(argv[0],"integrationPoints") == 0)
02112     theResponse = new ElementResponse(this, 10, Vector(numSections));
02113 
02114   else if (strcmp(argv[0],"integrationWeights") == 0)
02115     theResponse = new ElementResponse(this, 11, Vector(numSections));
02116 
02117   // section response -
02118   else if (strstr(argv[0],"section") != 0) {
02119     if (argc > 2) {
02120       int sectionNum = atoi(argv[1]);
02121       if (sectionNum > 0 && sectionNum <= numSections) {
02122 
02123         double xi[maxNumSections];
02124         double L = crdTransf->getInitialLength();
02125         beamIntegr->getSectionLocations(numSections, L, xi);
02126 
02127         output.tag("GaussPointOutput");
02128         output.attr("number",sectionNum);
02129         output.attr("eta",xi[sectionNum-1]*L);
02130 
02131         theResponse = sections[sectionNum-1]->setResponse(&argv[2], argc-2, eleInformation, output);
02132         
02133         output.endTag();
02134       }
02135     }
02136   }
02137   
02138   output.endTag(); // ElementOutput
02139 
02140   return theResponse;
02141 }
02142 
02143 int 
02144 ForceBeamColumn2d::getResponse(int responseID, Information &eleInfo)
02145 {
02146   static Vector vp(3);
02147   static Matrix fe(3,3);
02148 
02149   if (responseID == 1)
02150     return eleInfo.setVector(this->getResistingForce());
02151   
02152   else if (responseID == 2) {
02153     double p0[3]; p0[0] = 0.0; p0[1] = 0.0; p0[2] = 0.0;
02154     if (numEleLoads > 0)
02155       this->computeReactions(p0);
02156     theVector(3) =  Se(0);
02157     theVector(0) = -Se(0)+p0[0];
02158     theVector(2) = Se(1);
02159     theVector(5) = Se(2);
02160     double V = (Se(1)+Se(2))/crdTransf->getInitialLength();
02161     theVector(1) =  V+p0[1];
02162     theVector(4) = -V+p0[2];
02163     return eleInfo.setVector(theVector);
02164   }
02165 
02166   // Chord rotation
02167   else if (responseID == 7) {
02168     return eleInfo.setVector(Se);
02169   }
02170       
02171   // Chord rotation
02172   else if (responseID == 3) {
02173     vp = crdTransf->getBasicTrialDisp();
02174     return eleInfo.setVector(vp);
02175   }
02176 
02177   // Plastic rotation
02178   else if (responseID == 4) {
02179     this->getInitialFlexibility(fe);
02180     vp = crdTransf->getBasicTrialDisp();
02181     vp.addMatrixVector(1.0, fe, Se, -1.0);
02182     return eleInfo.setVector(vp);
02183   }
02184 
02185   // Point of inflection
02186   else if (responseID == 5) {
02187     double LI = 0.0;
02188 
02189     if (fabs(Se(1)+Se(2)) > DBL_EPSILON) {
02190       double L = crdTransf->getInitialLength();
02191       
02192       LI = Se(1)/(Se(1)+Se(2))*L;
02193     }
02194 
02195     return eleInfo.setDouble(LI);
02196   }
02197 
02198   // Tangent drift
02199   else if (responseID == 6) {
02200     double d2 = 0.0;
02201     double d3 = 0.0;
02202 
02203     double L = crdTransf->getInitialLength();
02204     
02205     // Location of inflection point from node I
02206     double LI = 0.0;
02207     if (fabs(Se(1)+Se(2)) > DBL_EPSILON)
02208       LI = Se(1)/(Se(1)+Se(2))*L;
02209       
02210     double wts[maxNumSections];
02211     beamIntegr->getSectionWeights(numSections, L, wts);
02212     
02213     double pts[maxNumSections];
02214     beamIntegr->getSectionLocations(numSections, L, pts);
02215     
02216     int i;
02217     for (i = 0; i < numSections; i++) {
02218       double x = pts[i]*L;
02219       if (x > LI)
02220         continue;
02221       const ID &type = sections[i]->getType();
02222       int order = sections[i]->getOrder();
02223       double kappa = 0.0;
02224       for (int j = 0; j < order; j++)
02225         if (type(j) == SECTION_RESPONSE_MZ)
02226           kappa += vs[i](j);
02227       double b = -LI+x;
02228       d2 += (wts[i]*L)*kappa*b;
02229     }
02230     
02231     d2 += beamIntegr->getTangentDriftI(L, LI, Se(1), Se(2));
02232     
02233     for (i = numSections-1; i >= 0; i--) {
02234       double x = pts[i]*L;
02235       if (x < LI)
02236         continue;
02237       const ID &type = sections[i]->getType();
02238       int order = sections[i]->getOrder();
02239       double kappa = 0.0;
02240       for (int j = 0; j < order; j++)
02241         if (type(j) == SECTION_RESPONSE_MZ)
02242           kappa += vs[i](j);
02243       double b = x-LI;
02244       d3 += (wts[i]*L)*kappa*b;
02245     }
02246     
02247     d3 += beamIntegr->getTangentDriftJ(L, LI, Se(1), Se(2));
02248 
02249     static Vector d(2);
02250     d(0) = d2;
02251     d(1) = d3;
02252 
02253     return eleInfo.setVector(d);
02254   }
02255 
02256   else if (responseID == 7)
02257     return eleInfo.setVector(Se);
02258 
02259   /*
02260   // Curvature sensitivity
02261   else if (responseID == 7) {
02262     Vector curv(numSections);
02263     for (int i = 0; i < numSections; i++) {
02264       int order = sections[i]->getOrder();
02265       const ID &type = sections[i]->getType();
02266       const Vector &dedh = sections[i]->getdedh();
02267       for (int j = 0; j < order; j++) {
02268         if (type(j) == SECTION_RESPONSE_MZ)
02269           curv(i) = dedh(j);
02270       }
02271     }
02272     return eleInfo.setVector(curv);
02273   }
02274 
02275   // Basic deformation sensitivity
02276   else if (responseID == 8) {  
02277     const Vector &dvdh = crdTransf->getBasicDisplSensitivity(1); // MHS hack
02278     return eleInfo.setVector(dvdh);
02279   }
02280   */
02281 
02282   // Plastic deformation sensitivity
02283   else if (responseID == 9) {  
02284     static Vector dvpdh(3);
02285 
02286     const Vector &dvdh = crdTransf->getBasicDisplSensitivity(1); // MHS hack
02287 
02288     dvpdh = dvdh;
02289     //opserr << dvpdh;
02290 
02291     static Matrix fe(3,3);
02292     this->getInitialFlexibility(fe);
02293 
02294     const Vector &dqdh = this->computedqdh(1); // MHS hack
02295 
02296     dvpdh.addMatrixVector(1.0, fe, dqdh, -1.0);
02297     //opserr << dvpdh;
02298 
02299     static Matrix fek(3,3);
02300     fek.addMatrixProduct(0.0, fe, kv, 1.0);
02301 
02302     dvpdh.addMatrixVector(1.0, fek, dvdh, -1.0);
02303     //opserr << dvpdh;
02304 
02305     const Matrix &dfedh = this->computedfedh(1); // MHS hack
02306 
02307     dvpdh.addMatrixVector(1.0, dfedh, Se, -1.0);
02308     //opserr << dvpdh << endln;
02309     //opserr << dfedh << endln;
02310 
02311     //opserr << dqdh + kv*dvdh << endln;
02312 
02313     return eleInfo.setVector(dvpdh);
02314   }
02315 
02316   else if (responseID == 10) {
02317     double L = crdTransf->getInitialLength();
02318     double pts[maxNumSections];
02319     beamIntegr->getSectionLocations(numSections, L, pts);
02320     Vector locs(numSections);
02321     for (int i = 0; i < numSections; i++)
02322       locs(i) = pts[i]*L;
02323     return eleInfo.setVector(locs);
02324   }
02325 
02326   else if (responseID == 11) {
02327     double L = crdTransf->getInitialLength();
02328     double wts[maxNumSections];
02329     beamIntegr->getSectionWeights(numSections, L, wts);
02330     Vector weights(numSections);
02331     for (int i = 0; i < numSections; i++)
02332       weights(i) = wts[i]*L;
02333     return eleInfo.setVector(weights);
02334   }
02335 
02336   else
02337     return -1;
02338 }
02339 
02340 int 
02341 ForceBeamColumn2d::getResponseSensitivity(int responseID, int gradNumber,
02342                                           Information &eleInfo)
02343 {
02344   // Basic deformation sensitivity
02345   if (responseID == 3) {  
02346     const Vector &dvdh = crdTransf->getBasicDisplSensitivity(gradNumber);
02347     return eleInfo.setVector(dvdh);
02348   }
02349 
02350   else if (responseID == 7) {
02351     static Vector dqdh(3);
02352 
02353     const Vector &dvdh = crdTransf->getBasicDisplSensitivity(gradNumber);
02354 
02355     dqdh.addMatrixVector(0.0, kv, dvdh, 1.0);
02356 
02357     dqdh.addVector(1.0, this->computedqdh(gradNumber), 1.0);
02358 
02359     return eleInfo.setVector(dqdh);
02360   }
02361 
02362   // Plastic deformation sensitivity
02363   else if (responseID == 4) {  
02364     static Vector dvpdh(3);
02365 
02366     const Vector &dvdh = crdTransf->getBasicDisplSensitivity(gradNumber);
02367 
02368     dvpdh = dvdh;
02369     //opserr << dvpdh;
02370 
02371     static Matrix fe(3,3);
02372     this->getInitialFlexibility(fe);
02373 
02374     const Vector &dqdh = this->computedqdh(gradNumber);
02375 
02376     dvpdh.addMatrixVector(1.0, fe, dqdh, -1.0);
02377     //opserr << dvpdh;
02378 
02379     static Matrix fek(3,3);
02380     fek.addMatrixProduct(0.0, fe, kv, 1.0);
02381 
02382     dvpdh.addMatrixVector(1.0, fek, dvdh, -1.0);
02383     //opserr << dvpdh;
02384 
02385     const Matrix &dfedh = this->computedfedh(gradNumber);
02386 
02387     dvpdh.addMatrixVector(1.0, dfedh, Se, -1.0);
02388     //opserr << dvpdh << endln;
02389     //opserr << dfedh << endln;
02390 
02391     //opserr << dqdh + kv*dvdh << endln;
02392 
02393     return eleInfo.setVector(dvpdh);
02394   }
02395 
02396   else
02397     return -1;
02398 }
02399 
02400 int
02401 ForceBeamColumn2d::setParameter(const char **argv, int argc, Parameter &param)
02402 {
02403   if (argc < 1)
02404     return -1;
02405 
02406   if (strcmp(argv[0],"rho") == 0)
02407     return param.addObject(1, this);
02408   
02409   // If the parameter belongs to a section or lower
02410   else if (strstr(argv[0],"section") != 0) {
02411     
02412     if (argc < 3)
02413       return -1;
02414 
02415     // Get section and material tag numbers from user input
02416     int paramSectionTag = atoi(argv[1]);
02417     
02418     // Find the right section and call its setParameter method
02419     int ok = 0;
02420     for (int i = 0; i < numSections; i++)
02421       if (paramSectionTag == sections[i]->getTag())
02422         ok += sections[i]->setParameter(&argv[2], argc-2, param);
02423 
02424     return ok;
02425   }
02426   
02427   else if (strstr(argv[0],"integration") != 0) {
02428     
02429     if (argc < 2)
02430       return -1;
02431 
02432     return beamIntegr->setParameter(&argv[1], argc-1, param);
02433   }
02434 
02435   else {
02436     opserr << "ForceBeamColumn2d::setParameter() - could not set parameter. " << endln;
02437     return -1;
02438   }
02439 }
02440 
02441 int
02442 ForceBeamColumn2d::updateParameter (int parameterID, Information &info)
02443 {
02444   if (parameterID == 1) {
02445     rho = info.theDouble;
02446     return 0;
02447   }
02448   else
02449     return -1;
02450 }
02451 
02452 int
02453 ForceBeamColumn2d::activateParameter(int passedParameterID)
02454 {
02455   parameterID = passedParameterID;
02456 
02457   return 0;  
02458 }
02459 
02460 const Matrix&
02461 ForceBeamColumn2d::getKiSensitivity(int gradNumber)
02462 {
02463   theMatrix.Zero();
02464   return theMatrix;
02465 }
02466 
02467 const Matrix&
02468 ForceBeamColumn2d::getMassSensitivity(int gradNumber)
02469 {
02470   theMatrix.Zero();
02471   return theMatrix;
02472 }
02473 
02474 const Vector&
02475 ForceBeamColumn2d::getResistingForceSensitivity(int gradNumber)
02476 {
02477   static Vector dqdh(3);
02478   dqdh = this->computedqdh(gradNumber);
02479 
02480   // Transform forces
02481   double dp0dh[3]; dp0dh[0] = 0.0; dp0dh[1] = 0.0; dp0dh[2] = 0.0;
02482   this->computeReactionSensitivity(dp0dh, gradNumber);
02483   Vector dp0dhVec(dp0dh, 3);
02484 
02485   static Vector P(6);
02486   P.Zero();
02487 
02488   if (crdTransf->isShapeSensitivity()) {
02489   // dAdh^T q
02490     P = crdTransf->getGlobalResistingForceShapeSensitivity(Se, dp0dhVec, gradNumber);
02491     // k dAdh u
02492     const Vector &dAdh_u = crdTransf->getBasicTrialDispShapeSensitivity();
02493     dqdh.addMatrixVector(1.0, kv, dAdh_u, 1.0);
02494   }
02495 
02496   // A^T (dqdh + k dAdh u)
02497   P += crdTransf->getGlobalResistingForce(dqdh, dp0dhVec);
02498 
02499   return P;
02500 }
02501 
02502 int
02503 ForceBeamColumn2d::commitSensitivity(int gradNumber, int numGrads)
02504 {
02505   int err = 0;
02506 
02507   double L = crdTransf->getInitialLength();
02508   double oneOverL = 1.0/L;
02509   
02510   double pts[maxNumSections];
02511   beamIntegr->getSectionLocations(numSections, L, pts);
02512   
02513   double wts[maxNumSections];
02514   beamIntegr->getSectionWeights(numSections, L, wts);
02515 
02516   double dLdh = crdTransf->getdLdh();
02517 
02518   double dptsdh[maxNumSections];
02519   beamIntegr->getLocationsDeriv(numSections, L, dLdh, dptsdh);
02520 
02521   double d1oLdh = crdTransf->getd1overLdh();
02522 
02523   static Vector dqdh(3);
02524   dqdh = this->computedqdh(gradNumber);
02525 
02526   // dvdh = A dudh + dAdh u
02527   const Vector &dvdh = crdTransf->getBasicDisplSensitivity(gradNumber);
02528   dqdh.addMatrixVector(1.0, kv, dvdh, 1.0);  // A dudh
02529 
02530   if (crdTransf->isShapeSensitivity()) {
02531     //const Vector &dAdh_u = crdTransf->getBasicTrialDispShapeSensitivity();
02532     //dqdh.addMatrixVector(1.0, kv, dAdh_u, 1.0);  // dAdh u
02533   }
02534 
02535   // Loop over integration points
02536   for (int i = 0; i < numSections; i++) {
02537 
02538     int order = sections[i]->getOrder();
02539     const ID &code = sections[i]->getType();
02540     
02541     double xL  = pts[i];
02542     double xL1 = xL-1.0;
02543 
02544     double dxLdh  = dptsdh[i];    
02545 
02546     Vector ds(workArea, order);
02547     int j;
02548     for (j = 0; j < order; j++) {
02549       switch(code(j)) {
02550       case SECTION_RESPONSE_P:
02551         ds(j) = dqdh(0);
02552         break;
02553       case SECTION_RESPONSE_MZ:
02554         ds(j) = xL1*dqdh(1) + xL*dqdh(2);
02555         break;
02556       case SECTION_RESPONSE_VY:
02557         ds(j) = oneOverL*(dqdh(1)+dqdh(2));
02558         break;
02559       default:
02560         ds(j) = 0.0;
02561         break;
02562       }
02563     }
02564 
02565     const Vector &dsdh = sections[i]->getStressResultantSensitivity(gradNumber,true);
02566     ds -= dsdh;
02567 
02568     for (j = 0; j < order; j++) {
02569       switch (code(j)) {
02570       case SECTION_RESPONSE_MZ:
02571         ds(j) += dxLdh*(Se(1)+Se(2));
02572         break;
02573       case SECTION_RESPONSE_VY:
02574         ds(j) += d1oLdh*(Se(1)+Se(2));
02575         break;
02576       default:
02577         break;
02578       }
02579     }
02580 
02581     Vector de(&workArea[order], order);
02582     const Matrix &fs = sections[i]->getSectionFlexibility();
02583     de.addMatrixVector(0.0, fs, ds, 1.0);
02584 
02585     err += sections[i]->commitSensitivity(de, gradNumber, numGrads);
02586   }
02587 
02588   return err;
02589 }
02590 
02591 const Vector &
02592 ForceBeamColumn2d::computedqdh(int gradNumber)
02593 {
02594   //opserr << "FBC2d::computedqdh " << gradNumber << endln;
02595 
02596   double L = crdTransf->getInitialLength();
02597   double oneOverL = 1.0/L;
02598   
02599   double pts[maxNumSections];
02600   beamIntegr->getSectionLocations(numSections, L, pts);
02601   
02602   double wts[maxNumSections];
02603   beamIntegr->getSectionWeights(numSections, L, wts);
02604 
02605   double dLdh = crdTransf->getdLdh();
02606 
02607   double dptsdh[maxNumSections];
02608   beamIntegr->getLocationsDeriv(numSections, L, dLdh, dptsdh);
02609 
02610   double dwtsdh[maxNumSections];
02611   beamIntegr->getWeightsDeriv(numSections, L, dLdh, dwtsdh);
02612 
02613   double d1oLdh = crdTransf->getd1overLdh();
02614 
02615   static Vector dvdh(3);
02616   dvdh.Zero();
02617 
02618   // Loop over the integration points
02619   for (int i = 0; i < numSections; i++) {
02620 
02621     int order = sections[i]->getOrder();
02622     const ID &code = sections[i]->getType();
02623     
02624     double xL  = pts[i];
02625     double xL1 = xL-1.0;
02626     double wtL = wts[i]*L;
02627     
02628     double dxLdh  = dptsdh[i];
02629     double dwtLdh = wts[i]*dLdh + dwtsdh[i]*L;
02630 
02631     //opserr << dptsdh[i] << ' ' << dwtsdh[i] << endln;
02632 
02633     // Get section stress resultant gradient
02634     Vector dsdh(&workArea[order], order);
02635     dsdh = sections[i]->getStressResultantSensitivity(gradNumber,true);
02636     //opserr << "FBC2d::dqdh -- " << gradNumber << ' ' << dsdh;
02637 
02638     // Add sensitivity wrt element loads
02639     if (numEleLoads > 0)
02640       this->computeSectionForceSensitivity(dsdh, i, gradNumber);
02641     
02642     int j;
02643     for (j = 0; j < order; j++) {
02644       switch (code(j)) {
02645       case SECTION_RESPONSE_MZ:
02646         dsdh(j) -= dxLdh*(Se(1)+Se(2));
02647         break;
02648       case SECTION_RESPONSE_VY:
02649         dsdh(j) -= d1oLdh*(Se(1)+Se(2));
02650         break;
02651       default:
02652         break;
02653       }
02654     }
02655 
02656     Vector dedh(workArea, order);
02657     const Matrix &fs = sections[i]->getSectionFlexibility();
02658     dedh.addMatrixVector(0.0, fs, dsdh, 1.0);
02659 
02660     for (j = 0; j < order; j++) {
02661       double dei = dedh(j)*wtL;
02662       switch(code(j)) {
02663       case SECTION_RESPONSE_P:
02664         dvdh(0) += dei; 
02665         break;
02666       case SECTION_RESPONSE_MZ:
02667         dvdh(1) += xL1*dei; 
02668         dvdh(2) += xL*dei;
02669         break;
02670       case SECTION_RESPONSE_VY:
02671         dei = oneOverL*dei;
02672         dvdh(1) += dei;
02673         dvdh(2) += dei;
02674       default:
02675         break;
02676       }
02677     }
02678 
02679     const Vector &e = vs[i];
02680     for (j = 0; j < order; j++) {
02681       switch(code(j)) {
02682       case SECTION_RESPONSE_P:
02683         dvdh(0) -= e(j)*dwtLdh;
02684         break;
02685       case SECTION_RESPONSE_MZ:
02686         dvdh(1) -= xL1*e(j)*dwtLdh;
02687         dvdh(2) -= xL*e(j)*dwtLdh;
02688         
02689         dvdh(1) -= dxLdh*e(j)*wtL;
02690         dvdh(2) -= dxLdh*e(j)*wtL;
02691         break;
02692       case SECTION_RESPONSE_VY:
02693         dvdh(1) -= oneOverL*e(j)*dwtLdh;
02694         dvdh(2) -= oneOverL*e(j)*dwtLdh;
02695 
02696         dvdh(1) -= d1oLdh*e(j)*wtL;
02697         dvdh(2) -= d1oLdh*e(j)*wtL;
02698         break;
02699       default:
02700         break;
02701       }
02702     }
02703   }
02704 
02705   static Matrix dfedh(3,3);
02706   dfedh.Zero();
02707 
02708   if (beamIntegr->addElasticFlexDeriv(L, dfedh, dLdh) < 0)
02709     dvdh.addMatrixVector(1.0, dfedh, Se, -1.0);
02710   
02711   //opserr << "dfedh: " << dfedh << endln;
02712 
02713   static Vector dqdh(3);
02714   dqdh.addMatrixVector(0.0, kv, dvdh, 1.0);
02715   
02716   return dqdh;
02717 }
02718 
02719 const Matrix&
02720 ForceBeamColumn2d::computedfedh(int gradNumber)
02721 {
02722   static Matrix dfedh(3,3);
02723 
02724   dfedh.Zero();
02725 
02726   double L = crdTransf->getInitialLength();
02727   double oneOverL  = 1.0/L;  
02728 
02729   double dLdh = crdTransf->getdLdh();
02730   double d1oLdh = crdTransf->getd1overLdh();
02731 
02732   beamIntegr->addElasticFlexDeriv(L, dfedh, dLdh);
02733 
02734   double xi[maxNumSections];
02735   beamIntegr->getSectionLocations(numSections, L, xi);
02736   
02737   double wt[maxNumSections];
02738   beamIntegr->getSectionWeights(numSections, L, wt);
02739 
02740   double dptsdh[maxNumSections];
02741   beamIntegr->getLocationsDeriv(numSections, L, dLdh, dptsdh);
02742 
02743   double dwtsdh[maxNumSections];
02744   beamIntegr->getWeightsDeriv(numSections, L, dLdh, dwtsdh);
02745 
02746   for (int i = 0; i < numSections; i++) {
02747 
02748     int order      = sections[i]->getOrder();
02749     const ID &code = sections[i]->getType();
02750     
02751     Matrix fb(workArea, order, NEBD);
02752     Matrix fb2(&workArea[order*NEBD], order, NEBD);
02753 
02754     double xL  = xi[i];
02755     double xL1 = xL-1.0;
02756     double wtL = wt[i]*L;
02757 
02758     double dxLdh  = dptsdh[i];
02759     double dwtLdh = wt[i]*dLdh + dwtsdh[i]*L;
02760 
02761     const Matrix &fs = sections[i]->getInitialFlexibility();
02762     const Matrix &dfsdh = sections[i]->getInitialFlexibilitySensitivity(gradNumber);
02763     fb.Zero();
02764     fb2.Zero();
02765 
02766     double tmp;
02767     int ii, jj;
02768     for (ii = 0; ii < order; ii++) {
02769       switch(code(ii)) {
02770       case SECTION_RESPONSE_P:
02771         for (jj = 0; jj < order; jj++) {
02772           fb(jj,0) += dfsdh(jj,ii)*wtL; // 1
02773 
02774           //fb(jj,0) += fs(jj,ii)*dwtLdh; // 3
02775 
02776           //fb2(jj,0) += fs(jj,ii)*wtL; // 4
02777         }
02778         break;
02779       case SECTION_RESPONSE_MZ:
02780         for (jj = 0; jj < order; jj++) {
02781           tmp = dfsdh(jj,ii)*wtL; // 1
02782           fb(jj,1) += xL1*tmp;
02783           fb(jj,2) += xL*tmp;
02784 
02785           tmp = fs(jj,ii)*wtL; // 2
02786           //fb(jj,1) += dxLdh*tmp;
02787           //fb(jj,2) += dxLdh*tmp;
02788 
02789           tmp = fs(jj,ii)*dwtLdh; // 3
02790           //fb(jj,1) += xL1*tmp;
02791           //fb(jj,2) += xL*tmp;
02792 
02793           tmp = fs(jj,ii)*wtL; // 4
02794           //fb2(jj,1) += xL1*tmp;
02795           //fb2(jj,2) += xL*tmp;
02796         }
02797         break;
02798       case SECTION_RESPONSE_VY:
02799         for (jj = 0; jj < order; jj++) {
02800           tmp = oneOverL*dfsdh(jj,ii)*wtL;
02801           fb(jj,1) += tmp;
02802           fb(jj,2) += tmp;
02803 
02804           // Need to complete for dLdh != 0
02805         }
02806         break;
02807       default:
02808         break;
02809       }
02810     }
02811     for (ii = 0; ii < order; ii++) {
02812       switch (code(ii)) {
02813       case SECTION_RESPONSE_P:
02814         for (jj = 0; jj < NEBD; jj++)
02815           dfedh(0,jj) += fb(ii,jj);
02816         break;
02817       case SECTION_RESPONSE_MZ:
02818         for (jj = 0; jj < NEBD; jj++) {
02819           tmp = fb(ii,jj); // 1,2,3
02820           dfedh(1,jj) += xL1*tmp;
02821           dfedh(2,jj) += xL*tmp;
02822 
02823           tmp = fb2(ii,jj); // 4
02824           //dfedh(1,jj) += dxLdh*tmp;
02825           //dfedh(2,jj) += dxLdh*tmp;
02826         }
02827         break;
02828       case SECTION_RESPONSE_VY:
02829         for (jj = 0; jj < NEBD; jj++) {
02830           tmp = oneOverL*fb(ii,jj);
02831           dfedh(1,jj) += tmp;
02832           dfedh(2,jj) += tmp;
02833 
02834           // Need to complete for dLdh != 0
02835         }
02836         break;
02837       default:
02838         break;
02839       }
02840     }
02841   }
02842   
02843   return dfedh;
02844 }
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