ForceBeamColumn3d.cpp

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

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