00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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
00059
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
00097
00098
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
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
00159
00160
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
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
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
00252 this->DomainComponent::setDomain(theDomain);
00253
00254
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
00264 if (crdTransf->initialize(theNodes[0], theNodes[1])) {
00265 opserr << "ForceBeamColumn3d::setDomain(): Error initializing coordinate transformation";
00266 exit(0);
00267 }
00268
00269
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
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
00301 if ((err = crdTransf->commitState()) != 0)
00302 return err;
00303
00304
00305 kvcommit = kv;
00306 Secommit = Se;
00307
00308
00309
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
00335 if ((err = crdTransf->revertToLastCommit()) != 0)
00336 return err;
00337
00338
00339 Se = Secommit;
00340 kv = kvcommit;
00341
00342 initialFlag = 0;
00343
00344
00345 return err;
00346 }
00347
00348 int ForceBeamColumn3d::revertToStart()
00349 {
00350
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
00366 if ((err = crdTransf->revertToStart()) != 0)
00367 return err;
00368
00369
00370 Secommit.Zero();
00371 kvcommit.Zero();
00372
00373 Se.Zero();
00374 kv.Zero();
00375
00376 initialFlag = 0;
00377
00378 return err;
00379 }
00380
00381
00382 const Matrix &
00383 ForceBeamColumn3d::getInitialStiff(void)
00384 {
00385
00386 if (Ki != 0)
00387 return *Ki;
00388
00389 static Matrix f(NEBD,NEBD);
00390 this->getInitialFlexibility(f);
00391
00392 static Matrix I(NEBD,NEBD);
00393 I.Zero();
00394 for (int i=0; i<NEBD; i++)
00395 I(i,i) = 1.0;
00396
00397
00398
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();
00412 return crdTransf->getGlobalStiffMatrix(kv, Se);
00413 }
00414
00415 const Vector &
00416 ForceBeamColumn3d::getResistingForce(void)
00417 {
00418
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
00441
00442 int
00443 ForceBeamColumn3d::update()
00444 {
00445
00446
00447 if (initialFlag == 2)
00448 this->revertToLastCommit();
00449
00450
00451 crdTransf->update();
00452
00453
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);
00475 static Matrix f(NEBD,NEBD);
00476
00477 static Matrix I(NEBD,NEBD);
00478 double dW;
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
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511 while (converged == false && numSubdivide <= maxSubdivisions) {
00512
00513
00514
00515
00516
00517 for (int l=0; l<3; l++) {
00518
00519
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
00529
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;
00538
00539 for (j=0; j <numIters; j++) {
00540
00541
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
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
00581
00582
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
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
00637 dSs = Ss;
00638 dSs.addVector(1.0, SsrSubdivide[i], -1.0);
00639
00640
00641 if (l == 0) {
00642
00643
00644
00645
00646 dvs.addMatrixVector(0.0, fsSubdivide[i], dSs, 1.0);
00647
00648 } else if (l == 2) {
00649
00650
00651
00652
00653
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
00665
00666
00667 const Matrix &fs0 = sections[i]->getInitialFlexibility();
00668 dvs.addMatrixVector(0.0, fs0, dSs, 1.0);
00669 }
00670
00671
00672 if (initialFlag != 0)
00673 vsSubdivide[i] += dvs;
00674
00675 sections[i]->setTrialSectionDeformation(vsSubdivide[i]);
00676
00677
00678 SsrSubdivide[i] = sections[i]->getStressResultant();
00679
00680
00681
00682 fsSubdivide[i] = sections[i]->getSectionFlexibility();
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695 dSs = Ss;
00696 dSs.addVector(1.0, SsrSubdivide[i], -1.0);
00697
00698 dvs.addMatrixVector(0.0, fsSubdivide[i], dSs, 1.0);
00699
00700
00701
00702
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
00794
00795
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
00833
00834
00835
00836 if (f.Solve(I, kvTrial) < 0)
00837 opserr << "ForceBeamColumn3d::update() -- could not invert flexibility\n";
00838
00839
00840
00841 dv = vin;
00842 dv += dvTrial;
00843 dv -= vr;
00844
00845
00846
00847
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
00857 if (fabs(dW) < tol) {
00858
00859
00860 dvToDo -= dvTrial;
00861 vin += dvTrial;
00862
00863
00864 if (dvToDo.Norm() <= DBL_EPSILON) {
00865 converged = true;
00866
00867 } else {
00868
00869
00870 dvTrial = dvToDo;
00871 numSubdivide = 1;
00872 }
00873
00874
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
00885 j = numIters+1;
00886 l = 4;
00887
00888 } else {
00889
00890
00891
00892 if (j == (numIters-1) && (l == 2)) {
00893 dvTrial /= factor;
00894 numSubdivide++;
00895 }
00896 }
00897
00898 }
00899 }
00900 }
00901 }
00902
00903
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
00912
00913
00914
00915
00916
00917
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:
00936 b(i,1) = xi - 1.0;
00937 b(i,2) = xi;
00938 break;
00939 case SECTION_RESPONSE_P:
00940 b(i,0) = 1.0;
00941 break;
00942 case SECTION_RESPONSE_VY:
00943 b(i,1) = b(i,2) = 1.0/L;
00944 break;
00945 case SECTION_RESPONSE_MY:
00946 b(i,3) = xi - 1.0;
00947 b(i,4) = xi;
00948 break;
00949 case SECTION_RESPONSE_VZ:
00950 b(i,3) = b(i,4) = 1.0/L;
00951 break;
00952 case SECTION_RESPONSE_T:
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:
00969 bp(i,1) = xi*(xi-1)*L*L/2;
00970 break;
00971 case SECTION_RESPONSE_P:
00972 bp(i,0) = (1-xi)*L;
00973 break;
00974 case SECTION_RESPONSE_VY:
00975 bp(i,1) = (xi-0.5)*L;
00976 break;
00977 case SECTION_RESPONSE_MY:
00978 bp(i,2) = xi*(1-xi)*L*L/2;
00979 break;
00980 case SECTION_RESPONSE_VZ:
00981 bp(i,2) = (0.5-xi)*L;
00982 break;
00983 case SECTION_RESPONSE_T:
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
01043 beamIntegr->addElasticDeformations(theLoad, loadFactor, L, v0);
01044
01045 if (type == LOAD_TAG_Beam3dUniformLoad) {
01046 double wy = data(0)*loadFactor;
01047 double wz = data(1)*loadFactor;
01048 double wx = data(2)*loadFactor;
01049
01050 Matrix &s_p = *sp;
01051
01052
01053 for (int i = 0; i < numSections; i++) {
01054 double x = xi[i]*L;
01055
01056 s_p(0,i) += wx*(L-x);
01057
01058 s_p(1,i) += wy*0.5*x*(x-L);
01059
01060 s_p(2,i) += wy*(x-0.5*L);
01061
01062 s_p(3,i) += wz*0.5*x*(L-x);
01063
01064 s_p(4,i) += wz*(x-0.5*L);
01065 }
01066
01067
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
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
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
01136 if (rho == 0.0)
01137 return 0;
01138
01139
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
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156 return 0;
01157 }
01158
01159 const Vector &
01160 ForceBeamColumn3d::getResistingForceIncInertia()
01161 {
01162
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
01180 if (alphaM != 0.0 || betaK != 0.0 || betaK0 != 0.0 || betaKc != 0.0)
01181 theVector += this->getRayleighDampingForces();
01182
01183 } else {
01184
01185
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
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
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
01248
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
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
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
01294 dData(loc++) = rho;
01295 dData(loc++) = tol;
01296
01297
01298
01299
01300
01301
01302 for (i=0; i<NEBD; i++)
01303 dData(loc++) = Secommit(i);
01304
01305
01306 for (i=0; i<NEBD; i++)
01307 for (j=0; j<NEBD; j++)
01308 dData(loc++) = kvcommit(i,j);
01309
01310
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
01329
01330 int dbTag = this->getDbTag();
01331 int i,j,k;
01332
01333 static ID idData(11);
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
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
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
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
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
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
01415
01416 if (numSections != idData(3)) {
01417
01418
01419
01420
01421
01422
01423
01424 if (numSections != 0) {
01425 for (int i=0; i<numSections; i++)
01426 delete sections[i];
01427 delete [] sections;
01428 }
01429
01430
01431 numSections = idData(3);
01432
01433
01434 if (vscommit != 0)
01435 delete [] vscommit;
01436
01437
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
01445 if (fs != 0)
01446 delete [] fs;
01447
01448
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
01456 if (vs != 0)
01457 delete [] vs;
01458
01459
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
01467 if (Ssr != 0)
01468 delete [] Ssr;
01469
01470
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
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
01512
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
01522 if (sections[i]->getClassTag() != sectClassTag) {
01523
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
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
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
01560 rho = dData(loc++);
01561 tol = dData(loc++);
01562
01563
01564
01565
01566
01567
01568 for (i=0; i<NEBD; i++)
01569 Secommit(i) = dData(loc++);
01570
01571
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
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
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
01722 static Vector ub(NEBD);
01723 ub = crdTransf->getBasicTrialDisp();
01724
01725 double L = crdTransf->getInitialLength();
01726
01727
01728 static double pts[maxNumSections];
01729 beamIntegr->getSectionLocations(numSections, L, pts);
01730
01731
01732 int i;
01733 double xi;
01734
01735
01736 Matrix ls(numSections, numSections);
01737 getCBDIinfluenceMatrix(numSections, pts, L, ls);
01738
01739
01740 Vector kappa_y(numSections);
01741 Vector kappa_z(numSections);
01742 static Vector vs;
01743
01744 for (i=0; i<numSections; i++) {
01745
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
01767 vs = sections[i]->getSectionDeformation();
01768
01769 kappa_z(i) = vs(sectionKey1);
01770 kappa_y(i) = vs(sectionKey2);
01771 }
01772
01773
01774
01775
01776 Vector v(numSections), w(numSections);
01777 static Vector xl(NDM), uxb(NDM);
01778 static Vector xg(NDM), uxg(NDM);
01779
01780
01781
01782 v.addMatrixVector (0.0, ls, kappa_z, 1.0);
01783
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
01795 sectionCoords[i] = crdTransf->getPointGlobalCoordFromLocal(xl);
01796
01797
01798
01799 uxb(0) = xi * ub(0);
01800 uxb(1) = v(i);
01801 uxb(2) = w(i);
01802
01803
01804 sectionDispls[i] = crdTransf->getPointGlobalDisplFromBasic(xi, uxb);
01805 }
01806 return;
01807 }
01808
01809 void
01810 ForceBeamColumn3d::Print(OPS_Stream &s, int flag)
01811 {
01812
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
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
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
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
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
01933 this->compSectionDisplacements(coords, displs);
01934
01935
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
01984
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
02033
02034
02035
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
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
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
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
02099 } else if (strcmp(argv[0],"inflectionPoint") == 0) {
02100 theResponse = new ElementResponse(this, 5, Vector(2));
02101
02102
02103 } else if (strcmp(argv[0],"tangentDrift") == 0) {
02104 theResponse = new ElementResponse(this, 6, Vector(4));
02105
02106
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
02141 double N = Se(0);
02142 theVector(6) = N;
02143 theVector(0) = -N+p0[0];
02144
02145
02146 double T = Se(5);
02147 theVector(9) = T;
02148 theVector(3) = -T;
02149
02150
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
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
02174 else if (responseID == 3) {
02175 vp = crdTransf->getBasicTrialDisp();
02176 return eleInfo.setVector(vp);
02177 }
02178
02179
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
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
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
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 ¶m)
02294 {
02295 if (argc < 1)
02296 return -1;
02297
02298
02299 if (strcmp(argv[0],"rho") == 0)
02300 return param.addObject(1, this);
02301
02302
02303 if (strstr(argv[0],"section") != 0) {
02304
02305 if (argc < 3)
02306 return -1;
02307
02308
02309 int paramSectionTag = atoi(argv[1]);
02310
02311
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
02337
02338
02339 return 0;
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368
02369
02370
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
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 }