BeamColumnJoint3d.cppGo 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.4 $ 00022 // $Date: 2006/08/04 22:22:37 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/element/joint/BeamColumnJoint3d.cpp,v $ 00024 00025 // Written: NM (nmitra@u.washington.edu) 00026 // Created: February 2003 00027 // Revised: September 2004 00028 // 00029 // Description: This file contains the class implementation for beam-column joint. 00030 // This element is a 4 noded 24 dof (6 dof at each node) finite area super-element, being a slight 00031 // variation of the 2d one. The element takes in 13 different material types in order to simulate 00032 // the inelastic action observed in a reinforced beam column joint. Though it has 6 dof per node 00033 // the out of the plane nodal dof are constrained or fixed and the inplane nodal dof are activated 00034 00035 #include <BeamColumnJoint3d.h> 00036 00037 #include <Domain.h> 00038 #include <Node.h> 00039 #include <Channel.h> 00040 #include <FEM_ObjectBroker.h> 00041 #include <Renderer.h> 00042 #include <MatrixUtil.h> 00043 #include <UniaxialMaterial.h> 00044 #include <string.h> 00045 #include <ElementResponse.h> 00046 00047 // class wide matrices 00048 Matrix BeamColumnJoint3d::Transf(12,24); 00049 Matrix BeamColumnJoint3d::Tran(3,6); 00050 00051 // full constructors: 00052 BeamColumnJoint3d::BeamColumnJoint3d(int tag,int Nd1, int Nd2, int Nd3, int Nd4, 00053 UniaxialMaterial& theMat1, 00054 UniaxialMaterial& theMat2, 00055 UniaxialMaterial& theMat3, 00056 UniaxialMaterial& theMat4, 00057 UniaxialMaterial& theMat5, 00058 UniaxialMaterial& theMat6, 00059 UniaxialMaterial& theMat7, 00060 UniaxialMaterial& theMat8, 00061 UniaxialMaterial& theMat9, 00062 UniaxialMaterial& theMat10, 00063 UniaxialMaterial& theMat11, 00064 UniaxialMaterial& theMat12, 00065 UniaxialMaterial& theMat13): 00066 Element(tag,ELE_TAG_BeamColumnJoint3d), 00067 connectedExternalNodes(4), elemActHeight(0.0), elemActWidth(0.0), 00068 elemWidth(0.0), elemHeight(0.0), HgtFac(1.0), WdtFac(1.0), 00069 Uecommit(24), UeIntcommit(4), UeprCommit(24), UeprIntCommit(4), 00070 BCJoint(13,16), dg_df(4,13), dDef_du(13,4), K(24,24), R(24), Node1(3), Node2(3), Node3(3), Node4(3) 00071 { 00072 // ensure the connectedExternalNode ID is of correct size & set values 00073 if (connectedExternalNodes.Size() != 4) 00074 opserr << "ERROR : BeamColumnJoint::BeamColumnJoint - " << tag << "failed to create an ID of size 4" << endln; 00075 00076 connectedExternalNodes(0) = Nd1 ; 00077 connectedExternalNodes(1) = Nd2 ; 00078 connectedExternalNodes(2) = Nd3 ; 00079 connectedExternalNodes(3) = Nd4 ; 00080 00081 MaterialPtr = new UniaxialMaterial*[13]; 00082 for (int x = 0; x <13; x++) 00083 { MaterialPtr[x] = 0; } 00084 00085 Uecommit.Zero(); 00086 UeIntcommit.Zero(); 00087 UeprCommit.Zero(); 00088 UeprIntCommit.Zero(); 00089 00090 BCJoint.Zero(); dg_df.Zero(); dDef_du.Zero(); 00091 K.Zero(); 00092 R.Zero(); 00093 Node1.Zero(); Node2.Zero(); Node3.Zero(); Node4.Zero(); 00094 00095 nodePtr[0] = 0; 00096 nodePtr[1] = 0; 00097 00098 // get a copy of the material and check we obtained a valid copy 00099 MaterialPtr[0] = theMat1.getCopy(); 00100 if (!MaterialPtr[0]){ 00101 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 1" << endln;} 00102 MaterialPtr[1] = theMat2.getCopy(); 00103 if (!MaterialPtr[1]){ 00104 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 2" << endln;} 00105 MaterialPtr[2] = theMat3.getCopy(); 00106 if (!MaterialPtr[2]){ 00107 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 3" << endln;} 00108 MaterialPtr[3] = theMat4.getCopy(); 00109 if (!MaterialPtr[3]){ 00110 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 4" << endln;} 00111 MaterialPtr[4] = theMat5.getCopy(); 00112 if (!MaterialPtr[4]){ 00113 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 5" << endln;} 00114 MaterialPtr[5] = theMat6.getCopy(); 00115 if (!MaterialPtr[5]){ 00116 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 6" << endln;} 00117 MaterialPtr[6] = theMat7.getCopy(); 00118 if (!MaterialPtr[6]){ 00119 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 7" << endln;} 00120 MaterialPtr[7] = theMat8.getCopy(); 00121 if (!MaterialPtr[7]){ 00122 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 8" << endln;} 00123 MaterialPtr[8] = theMat9.getCopy(); 00124 if (!MaterialPtr[8]){ 00125 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 9" << endln;} 00126 MaterialPtr[9] = theMat10.getCopy(); 00127 if (!MaterialPtr[9]){ 00128 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 10" << endln;} 00129 MaterialPtr[10] = theMat11.getCopy(); 00130 if (!MaterialPtr[10]){ 00131 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 11" << endln;} 00132 MaterialPtr[11] = theMat12.getCopy(); 00133 if (!MaterialPtr[11]){ 00134 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 12" << endln;} 00135 MaterialPtr[12] = theMat13.getCopy(); 00136 if (!MaterialPtr[12]){ 00137 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 13" << endln;} 00138 00139 } 00140 00141 // full constructors: 00142 BeamColumnJoint3d::BeamColumnJoint3d(int tag,int Nd1, int Nd2, int Nd3, int Nd4, 00143 UniaxialMaterial& theMat1, 00144 UniaxialMaterial& theMat2, 00145 UniaxialMaterial& theMat3, 00146 UniaxialMaterial& theMat4, 00147 UniaxialMaterial& theMat5, 00148 UniaxialMaterial& theMat6, 00149 UniaxialMaterial& theMat7, 00150 UniaxialMaterial& theMat8, 00151 UniaxialMaterial& theMat9, 00152 UniaxialMaterial& theMat10, 00153 UniaxialMaterial& theMat11, 00154 UniaxialMaterial& theMat12, 00155 UniaxialMaterial& theMat13, 00156 double elHgtFac, double elWdtFac): 00157 Element(tag,ELE_TAG_BeamColumnJoint3d), 00158 connectedExternalNodes(4), elemActHeight(0.0), elemActWidth(0.0), 00159 elemWidth(0), elemHeight(0), HgtFac(elHgtFac), WdtFac(elWdtFac), 00160 Uecommit(24), UeIntcommit(4), UeprCommit(24), UeprIntCommit(4), 00161 BCJoint(13,16), dg_df(4,13), dDef_du(13,4), K(24,24), R(24), Node1(3), Node2(3), Node3(3), Node4(3) 00162 { 00163 // ensure the connectedExternalNode ID is of correct size & set values 00164 if (connectedExternalNodes.Size() != 4) 00165 opserr << "ERROR : BeamColumnJoint::BeamColumnJoint - " << tag << "failed to create an ID of size 4" << endln; 00166 00167 connectedExternalNodes(0) = Nd1 ; 00168 connectedExternalNodes(1) = Nd2 ; 00169 connectedExternalNodes(2) = Nd3 ; 00170 connectedExternalNodes(3) = Nd4 ; 00171 00172 MaterialPtr = new UniaxialMaterial*[13]; 00173 for (int x = 0; x <13; x++) 00174 { MaterialPtr[x] = 0; } 00175 00176 Uecommit.Zero(); 00177 UeIntcommit.Zero(); 00178 UeprCommit.Zero(); 00179 UeprIntCommit.Zero(); 00180 00181 BCJoint.Zero(); dg_df.Zero(); dDef_du.Zero(); 00182 K.Zero(); 00183 R.Zero(); 00184 Node1.Zero(); Node2.Zero(); Node3.Zero(); Node4.Zero(); 00185 00186 nodePtr[0] = 0; 00187 nodePtr[1] = 0; 00188 00189 // get a copy of the material and check we obtained a valid copy 00190 MaterialPtr[0] = theMat1.getCopy(); 00191 if (!MaterialPtr[0]){ 00192 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 1" << endln;} 00193 MaterialPtr[1] = theMat2.getCopy(); 00194 if (!MaterialPtr[1]){ 00195 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 2" << endln;} 00196 MaterialPtr[2] = theMat3.getCopy(); 00197 if (!MaterialPtr[2]){ 00198 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 3" << endln;} 00199 MaterialPtr[3] = theMat4.getCopy(); 00200 if (!MaterialPtr[3]){ 00201 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 4" << endln;} 00202 MaterialPtr[4] = theMat5.getCopy(); 00203 if (!MaterialPtr[4]){ 00204 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 5" << endln;} 00205 MaterialPtr[5] = theMat6.getCopy(); 00206 if (!MaterialPtr[5]){ 00207 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 6" << endln;} 00208 MaterialPtr[6] = theMat7.getCopy(); 00209 if (!MaterialPtr[6]){ 00210 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 7" << endln;} 00211 MaterialPtr[7] = theMat8.getCopy(); 00212 if (!MaterialPtr[7]){ 00213 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 8" << endln;} 00214 MaterialPtr[8] = theMat9.getCopy(); 00215 if (!MaterialPtr[8]){ 00216 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 9" << endln;} 00217 MaterialPtr[9] = theMat10.getCopy(); 00218 if (!MaterialPtr[9]){ 00219 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 10" << endln;} 00220 MaterialPtr[10] = theMat11.getCopy(); 00221 if (!MaterialPtr[10]){ 00222 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 11" << endln;} 00223 MaterialPtr[11] = theMat12.getCopy(); 00224 if (!MaterialPtr[11]){ 00225 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 12" << endln;} 00226 MaterialPtr[12] = theMat13.getCopy(); 00227 if (!MaterialPtr[12]){ 00228 opserr << "ERROR : BeamColumnJoint::Constructor failed to get a copy of material 13" << endln;} 00229 00230 } 00231 00232 // default constructor: 00233 BeamColumnJoint3d::BeamColumnJoint3d(): 00234 Element(0,ELE_TAG_BeamColumnJoint3d), 00235 connectedExternalNodes(4), elemActHeight(0.0), elemActWidth(0.0), 00236 elemWidth(0), elemHeight(0), HgtFac(1.0), WdtFac(1.0), 00237 Uecommit(24), UeIntcommit(4), UeprCommit(24), UeprIntCommit(4), 00238 BCJoint(13,16), dg_df(4,13), dDef_du(13,4), K(24,24), R(24), Node1(3), Node2(3), Node3(3), Node4(3) 00239 { 00240 nodePtr[0] = 0; 00241 nodePtr[1] = 0; 00242 for (int x = 0; x <13; x++) 00243 { MaterialPtr[x] = 0; } 00244 // does nothing (invoked by FEM_ObjectBroker) 00245 } 00246 00247 // destructor: 00248 BeamColumnJoint3d::~BeamColumnJoint3d() 00249 { 00250 for (int i =0; i<13; i++) 00251 { 00252 if (MaterialPtr[i] != 0) 00253 delete MaterialPtr[i]; 00254 } 00255 00256 if (MaterialPtr) 00257 delete [] MaterialPtr; 00258 } 00259 00260 // public methods 00261 00262 int 00263 BeamColumnJoint3d::getNumExternalNodes(void) const 00264 { 00265 return 4; 00266 } 00267 00268 const ID & 00269 BeamColumnJoint3d::getExternalNodes(void) 00270 { 00271 return connectedExternalNodes; 00272 } 00273 00274 Node **BeamColumnJoint3d::getNodePtrs(void) 00275 { 00276 return nodePtr; 00277 } 00278 00279 int 00280 BeamColumnJoint3d::getNumDOF(void) 00281 { 00282 return 24; 00283 } 00284 00285 void 00286 BeamColumnJoint3d::setDomain(Domain *theDomain) 00287 { 00288 if (theDomain == 0) 00289 opserr << "ERROR : BeamColumnJoint::setDomain -- Domain is null" << endln; 00290 00291 // Check Domain is not null - invoked when object removed from a domain 00292 if (theDomain == 0) { 00293 nodePtr[0] = 0; 00294 nodePtr[1] = 0; 00295 } 00296 00297 //node pointers set 00298 for (int i = 0; i < 4; i++ ) { 00299 nodePtr[i] = theDomain->getNode( connectedExternalNodes(i) ) ; 00300 if (nodePtr[i] == 0) 00301 { 00302 opserr << "ERROR : BeamColumnJoint::setDomain -- node pointer is null" << endln; 00303 exit(-1); // donot go any further - otherwise segmentation fault 00304 } 00305 } 00306 00307 // call the base class method 00308 this->DomainComponent::setDomain(theDomain); 00309 00310 // ensure connected nodes have correct dof's 00311 int dofNd1 = nodePtr[0]->getNumberDOF(); 00312 int dofNd2 = nodePtr[1]->getNumberDOF(); 00313 int dofNd3 = nodePtr[2]->getNumberDOF(); 00314 int dofNd4 = nodePtr[3]->getNumberDOF(); 00315 00316 // check for proper nodal degrees of freedom 00317 if ((dofNd1 != 6) || (dofNd2 != 6) || (dofNd3 != 6) || (dofNd4 != 6)) 00318 { 00319 opserr << "ERROR : BeamColumnJoint::setDomain -- number of DOF associated with the node incorrect" << endln; 00320 exit(-1); // donot go any further - otherwise segmentation fault 00321 } 00322 00323 // obtain the nodal coordinates 00324 const Vector &end1Crd = nodePtr[0]->getCrds(); 00325 const Vector &end2Crd = nodePtr[1]->getCrds(); 00326 const Vector &end3Crd = nodePtr[2]->getCrds(); 00327 const Vector &end4Crd = nodePtr[3]->getCrds(); 00328 00329 // determine the element width and the element height. It should be noted that the element width is determined by nodal coordinates 2 and 4 00330 // whereas the element height is determined by the nodal coordinates 1 and 3 00331 // general order anticlockwise arrangement , but can be clockwise too. 00332 Node1 = end1Crd; 00333 Node2 = end2Crd; 00334 Node3 = end3Crd; 00335 Node4 = end4Crd; 00336 00337 Node3 = Node3 - Node1; 00338 Node2 = Node2 - Node4; 00339 00340 double beam = 1.0; //HgtFac; 00341 double col = 1.0; //WdtFac; 00342 00343 elemActHeight = fabs(Node3.Norm()); 00344 elemActWidth = fabs(Node2.Norm()); 00345 elemHeight = beam*elemActHeight; 00346 elemWidth = col*elemActWidth; 00347 00348 if ((elemHeight <= 1e-12) || (elemWidth <= 1e-12)) 00349 { 00350 opserr << "ERROR : BeamColumnJoint::setDomain -- length or width not correct, division by zero occurs"<< endln; 00351 exit(-1); // donot go any further - otherwise segmentation fault 00352 } 00353 00354 getBCJoint(); 00355 getdg_df(); 00356 getdDef_du(); 00357 formTransfMat(); 00358 } 00359 00360 int 00361 BeamColumnJoint3d::commitState() 00362 { 00363 // store commited external nodal displacements 00364 Uecommit = UeprCommit; 00365 // store commited internal nodal displacements 00366 UeIntcommit = UeprIntCommit; 00367 00368 // store material history data. 00369 int mcs = 0; 00370 for (int j=0; j<13; j++) 00371 { 00372 if (MaterialPtr[j] != 0) mcs = MaterialPtr[j]->commitState(); 00373 if (mcs != 0) break; 00374 } 00375 00376 return mcs; 00377 } 00378 00379 int 00380 BeamColumnJoint3d::revertToLastCommit() 00381 { 00382 00383 int mcs = 0; 00384 for (int j=0; j<13; j++) 00385 { 00386 if (MaterialPtr[j] != 0) mcs = MaterialPtr[j]->revertToLastCommit(); 00387 if (mcs != 0) break; 00388 } 00389 UeprCommit = Uecommit; 00390 UeprIntCommit = UeIntcommit; 00391 00392 this->update(); 00393 00394 return mcs; 00395 } 00396 00397 int 00398 BeamColumnJoint3d::revertToStart() 00399 { 00400 int mcs = 0; 00401 for (int j=0; j<13; j++) 00402 { 00403 if (MaterialPtr[j] != 0) mcs = MaterialPtr[j]->revertToStart(); 00404 if (mcs != 0) break; 00405 } 00406 00407 return mcs; 00408 } 00409 00410 int 00411 BeamColumnJoint3d::update() 00412 { 00413 00414 Vector Ue(28); 00415 Ue.Zero(); 00416 00417 // determine commited displacements given trial displacements 00418 this->getGlobalDispls(Ue); 00419 00420 // update displacements for the external nodes 00421 UeprCommit.Extract(Ue,0,1.0); 00422 00423 // update displacement for the internal nodes 00424 UeprIntCommit.Extract(Ue,24,1.0); 00425 00426 return 0; 00427 } 00428 00429 00430 const Matrix & 00431 BeamColumnJoint3d::getTangentStiff(void) 00432 { 00433 return K; 00434 } 00435 00436 const Matrix & 00437 BeamColumnJoint3d::getInitialStiff(void) 00438 { 00439 return getTangentStiff(); 00440 } 00441 00442 const Vector & 00443 BeamColumnJoint3d::getResistingForce() 00444 { 00445 return R; 00446 } 00447 00448 void BeamColumnJoint3d::getGlobalDispls(Vector &dg) 00449 { 00450 // local variables that will be used in this method 00451 int converge = 0; 00452 int linesearch = 0; 00453 int totalCount = 0; 00454 int dtConverge = 0; 00455 int incCount = 0; 00456 int count = 0; 00457 int maxTotalCount = 1000; 00458 int maxCount = 20; 00459 double loadStep = 0.0; 00460 double dLoadStep = 1.0; 00461 double stepSize; 00462 00463 Vector uExtOld(24); uExtOld.Zero(); 00464 Vector uExt(12); uExt.Zero(); 00465 Vector duExt(12); duExt.Zero(); 00466 Vector uIntOld(4); uIntOld.Zero(); 00467 Vector uInt(4); uInt.Zero(); 00468 Vector duInt(4); duInt.Zero(); 00469 Vector duIntTemp(4); duIntTemp.Zero(); 00470 Vector intEq(4); intEq.Zero(); 00471 Vector intEqLast(4); intEqLast.Zero(); 00472 Vector Uepr(24); Uepr.Zero(); 00473 Vector UeprInt(4); UeprInt.Zero(); 00474 Vector Ut(24); Ut.Zero(); 00475 Vector duExtTemp(24); duExtTemp.Zero(); 00476 00477 Vector disp1 = nodePtr[0]->getTrialDisp(); 00478 Vector disp2 = nodePtr[1]->getTrialDisp(); 00479 Vector disp3 = nodePtr[2]->getTrialDisp(); 00480 Vector disp4 = nodePtr[3]->getTrialDisp(); 00481 00482 for (int i = 0; i < 6; i++) 00483 { 00484 Ut(i) = disp1(i); 00485 Ut(i+6) = disp2(i); 00486 Ut(i+12) = disp3(i); 00487 Ut(i+18) = disp4(i); 00488 } 00489 00490 Uepr = Uecommit; 00491 00492 UeprInt = UeIntcommit; 00493 00494 uExtOld = Uepr; 00495 00496 duExtTemp = Ut - Uepr; 00497 duExt.addMatrixVector(0.0,Transf,duExtTemp,1.0); 00498 uExt.addMatrixVector(0.0,Transf,uExtOld,1.0); 00499 00500 uIntOld = UeprInt; 00501 uInt = uIntOld; 00502 00503 double tol = 1e-12; 00504 double tolIntEq = tol; 00505 double toluInt = (tol>tol*uInt.Norm())? tol:tol*uInt.Norm(); 00506 double tolIntEqdU = tol; 00507 double ctolIntEqdU = tol; 00508 double ctolIntEq = tol; 00509 double normDuInt = toluInt; 00510 double normIntEq = tolIntEq; 00511 double normIntEqdU = tolIntEqdU; 00512 00513 Vector u(16); u.Zero(); 00514 00515 double engrLast = 0.0; 00516 double engr = 0.0; 00517 00518 Vector fSpring(13); fSpring.Zero(); 00519 Vector kSpring(13); kSpring.Zero(); 00520 Matrix dintEq_du(4,4); dintEq_du.Zero(); 00521 Matrix df_dDef(13,13); df_dDef.Zero(); 00522 Matrix tempintEq_du (4,13); tempintEq_du.Zero(); 00523 00524 00525 while ((loadStep < 1.0) && (totalCount < maxTotalCount)) 00526 { 00527 count = 0; 00528 converge = 0; 00529 dtConverge = 0; 00530 while ((!converge) && (count < maxCount)) 00531 { 00532 if (dLoadStep <= 1e-3) 00533 { 00534 dLoadStep = dLoadStep; 00535 } 00536 totalCount ++; 00537 count ++; 00538 00539 for (int ic = 0; ic < 12; ic++ ) { 00540 u(ic) = uExt(ic) + duExt(ic); 00541 } 00542 u(12) = uInt(0); 00543 u(13) = uInt(1); 00544 u(14) = uInt(2); 00545 u(15) = uInt(3); 00546 00547 getMatResponse(u,fSpring,kSpring); 00548 00549 00550 // performs internal equilibrium 00551 00552 intEq(0) = -fSpring(2)-fSpring(3)+fSpring(9)-fSpring(12)/elemHeight; 00553 intEq(1) = fSpring(1)-fSpring(5)-fSpring(7)+fSpring(12)/elemWidth; 00554 intEq(2) = -fSpring(4)-fSpring(8)+fSpring(10)+fSpring(12)/elemHeight; 00555 intEq(3) = fSpring(0)-fSpring(6)-fSpring(11)-fSpring(12)/elemWidth; 00556 00557 matDiag(kSpring, df_dDef); 00558 00560 tempintEq_du.addMatrixProduct(0.0,dg_df,df_dDef,1.0); 00561 dintEq_du.addMatrixProduct(0.0,tempintEq_du,dDef_du,1.0); 00562 normIntEq = intEq.Norm(); 00563 normIntEqdU = 0.0; 00564 for (int jc = 0; jc<4 ; jc++) 00565 { 00566 normIntEqdU += intEq(jc)*duInt(jc); 00567 } 00568 normIntEqdU = fabs(normIntEqdU); 00569 00570 if (totalCount == 1) 00571 { 00572 tolIntEq = (tol>tol*normIntEq) ? tol:tol*normIntEq; 00573 tolIntEqdU = tol; 00574 } 00575 else if (totalCount == 2) 00576 { 00577 tolIntEqdU = (tol>tol*normIntEqdU) ? tol:tol*normIntEqdU; 00578 } 00579 ctolIntEqdU = (tolIntEqdU*dLoadStep > tol) ? tolIntEqdU*dLoadStep:tol; 00580 ctolIntEq = (tolIntEq*dLoadStep > tol) ? tolIntEq*dLoadStep:tol; 00581 00582 // check for convergence starts 00583 if ((normIntEq < tol) || ((normIntEqdU < tol) && (count >1)) || (normDuInt < toluInt) || (dLoadStep < 1e-3)) 00584 { 00585 if ((normIntEq > ctolIntEq) || (normIntEqdU > tolIntEqdU) || (normDuInt > toluInt)) 00586 { 00587 dtConverge = 1; 00588 } 00589 else 00590 { 00591 dtConverge = 0; 00592 } 00593 00594 converge = 1; 00595 loadStep = loadStep + dLoadStep; 00596 if (fabs(1.0 - loadStep) < tol) 00597 { 00598 loadStep = 1.0; 00599 } 00600 } 00601 else 00602 { 00604 dintEq_du.Solve(intEq,duInt); 00605 duInt *= -1; 00606 00607 normDuInt = duInt.Norm(); 00608 if (!linesearch) 00609 { 00610 uInt = uInt + duInt; 00611 } 00612 else 00613 { 00614 engrLast = 0.0; 00615 engr = 0.0; 00616 00617 for (int jd = 0; jd<4 ; jd++) 00618 { 00619 engrLast += duInt(jd)*intEqLast(jd); 00620 engr += duInt(jd)*intEq(jd); 00621 } 00622 00623 if (fabs(engr) > tol*engrLast) 00624 { 00625 duIntTemp = duInt; 00626 duIntTemp *= -1; 00627 // lineSearch algorithm requirement 00628 stepSize = getStepSize(engrLast,engr,uExt,duExt,uInt,duIntTemp,tol); 00629 00630 if (fabs(stepSize) > 0.001) 00631 { 00632 uInt = uInt + stepSize*duInt; 00633 } 00634 else 00635 { 00636 uInt = uInt + duInt; 00637 } 00638 } 00639 else 00640 { 00641 uInt = uInt + duInt; 00642 } 00643 intEqLast = intEq; 00644 } 00645 } 00646 } 00647 00648 if (!converge && loadStep < 1.0) 00649 { 00650 incCount = 0; 00651 maxCount = 25; 00652 if (!linesearch) 00653 { 00654 linesearch = 1; 00655 uInt = uIntOld; 00656 duInt.Zero(); 00657 } 00658 else 00659 { 00660 opserr << "WARNING : BeamColumnJoint::getGlobalDispls() - convergence problem in state determination" << endln; 00661 00662 uInt = uIntOld; 00663 duInt.Zero(); 00664 duExt = duExt*0.1; 00665 00666 dLoadStep = dLoadStep*0.1; 00667 } 00668 } 00669 else if (loadStep < 1.0) 00670 { 00671 maxCount = 10; 00672 incCount ++; 00673 normDuInt = toluInt; 00674 if ((incCount < maxCount) || dtConverge) 00675 { 00676 uExt = uExt + duExt; 00677 if (loadStep + dLoadStep > 1.0) 00678 { 00679 duExt = duExt*(1.0 - loadStep)/dLoadStep; 00680 dLoadStep = 1.0 - loadStep; 00681 incCount = 9; 00682 } 00683 } 00684 else 00685 { 00686 incCount = 0; 00687 uExt = uExt + duExt; 00688 dLoadStep = dLoadStep*10; 00689 if (loadStep + dLoadStep > 1.0) 00690 { 00691 uExt = uExt + duExt*(1.0 - loadStep)/dLoadStep; 00692 dLoadStep = 1.0 - loadStep; 00693 incCount = 9; 00694 } 00695 } 00696 } 00697 } 00698 00699 00700 // determination of stiffness matrix and the residual force vector for the element 00701 00702 formR(fSpring); 00703 00704 formK(kSpring); 00705 00706 for (int ig = 0; ig < 25; ig++ ) { 00707 if (ig<24) 00708 { 00709 dg(ig) = Ut(ig); 00710 } 00711 } 00712 00713 00714 dg(24) = uInt(0); 00715 dg(25) = uInt(1); 00716 dg(26) = uInt(2); 00717 dg(27) = uInt(3); 00718 00719 } 00720 00721 void BeamColumnJoint3d::getMatResponse(Vector U, Vector &fS, Vector &kS) 00722 { 00723 double jh = HgtFac; // factor for beams 00724 double jw = WdtFac; // factor for column 00725 00726 // obtains the material response from the material class 00727 Vector defSpring(13); 00728 defSpring.Zero(); 00729 fS.Zero(); 00730 kS.Zero(); 00731 00732 defSpring.addMatrixVector(0.0, BCJoint, U, 1.0); 00733 00734 // slip @ bar = slip @ spring * tc couple 00735 00736 defSpring(0) = defSpring(0)*jw; // new model applied on 00737 defSpring(1) = defSpring(1)*jw; // 4th June 2004 00738 defSpring(6) = defSpring(6)*jw; // h, h_bar distinction was 00739 defSpring(7) = defSpring(7)*jw; // removed from the previous model 00740 00741 defSpring(3) = defSpring(3)*jh; // by making h = h_bar 00742 defSpring(4) = defSpring(4)*jh; // similar case done for width 00743 defSpring(9) = defSpring(9)*jh; 00744 defSpring(10) = defSpring(10)*jh; 00745 00746 for (int j=0; j<13; j++) 00747 { 00748 MaterialPtr[j]->setTrialStrain(defSpring[j]); 00749 kS[j] = MaterialPtr[j]->getTangent(); 00750 fS[j] = MaterialPtr[j]->getStress(); 00751 } 00752 00753 // force @ spring = force @ bar * tc couple 00754 00755 fS(0) = fS(0)*jw; 00756 fS(1) = fS(1)*jw; 00757 fS(6) = fS(6)*jw; 00758 fS(7) = fS(7)*jw; 00759 00760 fS(3) = fS(3)*jh; 00761 fS(4) = fS(4)*jh; 00762 fS(9) = fS(9)*jh; 00763 fS(10) = fS(10)*jh; 00764 00765 // stiffness @ spring = stiffness @ bar * (tc couple)^2 00766 00767 kS(0) = kS(0)*jw*jw; 00768 kS(1) = kS(1)*jw*jw; 00769 kS(6) = kS(6)*jw*jw; 00770 kS(7) = kS(7)*jw*jw; 00771 00772 kS(3) = kS(3)*jh*jh; 00773 kS(4) = kS(4)*jh*jh; 00774 kS(9) = kS(9)*jh*jh; 00775 kS(10) = kS(10)*jh*jh; 00776 } 00777 00778 00779 void BeamColumnJoint3d::getdDef_du() 00780 { 00781 00782 dDef_du.Zero(); 00783 for (int jb=0; jb<13; jb++) 00784 { 00785 dDef_du(jb,0) = BCJoint(jb,12); 00786 dDef_du(jb,1) = BCJoint(jb,13); 00787 dDef_du(jb,2) = BCJoint(jb,14); 00788 dDef_du(jb,3) = BCJoint(jb,15); 00789 } 00790 } 00791 00792 void BeamColumnJoint3d::matDiag(Vector k,Matrix &dfd) 00793 { 00794 dfd.Zero(); 00795 // takes in a vector and converts it to a diagonal matrix (could have been placed as a method in matrix class) 00796 for (int ja=0; ja<13; ja++) 00797 { 00798 dfd(ja,ja) = k(ja); 00799 } 00800 00801 } 00802 00803 void BeamColumnJoint3d::formR(Vector f) 00804 { 00805 00806 // develops the element residual force vector 00807 Vector rForceTemp(16); 00808 Vector Rtempo(12); 00809 rForceTemp.Zero(); 00810 00811 rForceTemp.addMatrixTransposeVector(0.0,BCJoint,f,1.0); // rForceTemp = BCJoint'*f 00812 Rtempo.Extract(rForceTemp,0,1.0); 00813 00814 R.addMatrixTransposeVector(0.0,Transf,Rtempo,1.0); // R = Transf'*Rtempo 00815 } 00816 00817 void BeamColumnJoint3d::formK(Vector k) 00818 { 00819 // develops the element stiffness matrix 00820 Matrix kSprDiag(13,13); 00821 kSprDiag.Zero(); 00822 Matrix kRForce(16,16); 00823 kRForce.Zero(); 00824 Matrix kRFT1(4,12); 00825 kRFT1.Zero(); 00826 Matrix kRFT2(4,4); 00827 kRFT2.Zero(); 00828 Matrix kRFT3(12,4); 00829 kRFT3.Zero(); 00830 Matrix I(4,4); 00831 I.Zero(); 00832 Matrix kRSTinv(4,4); 00833 kRSTinv.Zero(); 00834 Matrix kRF(12,12); 00835 kRF.Zero(); 00836 Matrix K2Temp(12,4); 00837 K2Temp.Zero(); 00838 Matrix K2(12,12); 00839 K2.Zero(); 00840 00841 matDiag(k,kSprDiag); 00842 00843 kRForce.addMatrixTripleProduct(0.0,BCJoint,kSprDiag,1.0); // kRForce = BCJoint'*kSprDiag*BCJoint 00844 kRFT2.Extract(kRForce,12,12,1.0); 00845 kRFT1.Extract(kRForce,12,0,1.0); 00846 kRFT3.Extract(kRForce,0,12,1.0); 00847 kRF.Extract(kRForce,0,0,1.0); 00848 00849 for (int ic=0; ic<4; ic++) 00850 { 00851 I(ic,ic) = 1.0; 00852 } 00853 kRFT2.Solve(I,kRSTinv); 00854 00855 K2Temp.addMatrixProduct(0.0,kRFT3,kRSTinv,1.0); 00856 00857 // done for some idiotic reason 00858 for(int i = 0; i < 12; ++i) 00859 { 00860 for(int j = 0; j < 4; ++j) 00861 { 00862 if(fabs(K2Temp(i,j)) < 1e-15) 00863 K2Temp(i,j) = 0.; 00864 } 00865 } 00866 00867 K2.addMatrixProduct(0.0,K2Temp,kRFT1,1.0); // K2 = kRFT3*kRSTinv*kRFT1 00868 00869 // done for some idiotic reason 00870 for(int i1 = 0; i1 < 12; ++i1) 00871 { 00872 for(int j1 = 0; j1 < 12; ++j1) 00873 { 00874 if(fabs(K2(i1,j1)) < 1e-15) 00875 K2(i1,j1) = 0.; 00876 } 00877 } 00878 00879 kRF.addMatrix(1.0,K2,-1.0); // kRF = kRF - K2 00880 00881 K.addMatrixTripleProduct(0.0,Transf,kRF,1.0); // K = Transf'*kRF*Transf 00882 } 00883 00884 void BeamColumnJoint3d::formTransfMat() 00885 { 00886 Transf.Zero(); 00887 Tran.Zero(); 00888 00889 double Crd24 = fabs(sqrt((Node2(0)*Node2(0)) + (Node2(1)*Node2(1)) + (Node2(2)*Node2(2)))); 00890 double Crd31 = fabs(sqrt((Node3(0)*Node3(0)) + (Node3(1)*Node3(1)) + (Node3(2)*Node3(2)))); 00891 00892 double a24 = Node2(0)/Crd24; double b24 = Node2(1)/Crd24; double c24 = Node2(2)/Crd24; 00893 double a31 = Node3(0)/Crd31; double b31 = Node3(1)/Crd31; double c31 = Node3(2)/Crd31; 00894 00895 Tran(0,0) = a24; Tran(0,1) = b24; Tran(0,2) = c24; 00896 Tran(1,0) = a31; Tran(1,1) = b31; Tran(1,2) = c31; 00897 Tran(2,3) = b24*c31 - b31*c24; Tran(2,4) = -a24*c31 + c24*a31; Tran(2,5) = a24*b31 - b24*a31; 00898 00899 Transf.Assemble(Tran,0,0,1.0); 00900 Transf.Assemble(Tran,3,6,1.0); 00901 Transf.Assemble(Tran,6,12,1.0); 00902 Transf.Assemble(Tran,9,18,1.0); 00903 00904 } 00905 00906 void BeamColumnJoint3d::getdg_df() 00907 { 00908 dg_df.Zero(); // basically derivative of intEq 00909 dg_df(0,2) = -1; 00910 dg_df(0,3) = -1; 00911 dg_df(0,9) = 1; 00912 dg_df(0,12) = -1/elemHeight; 00913 dg_df(1,1) = 1; 00914 dg_df(1,5) = -1; 00915 dg_df(1,7) = -1; 00916 dg_df(1,12) = 1/elemWidth; 00917 dg_df(2,4) = -1; 00918 dg_df(2,8) = -1; 00919 dg_df(2,10) = 1; 00920 dg_df(2,12) = 1/elemHeight; 00921 dg_df(3,0) = 1; 00922 dg_df(3,6) = -1; 00923 dg_df(3,11) = -1; 00924 dg_df(3,12) = -1/elemWidth; 00925 00926 } 00927 00928 void BeamColumnJoint3d::getBCJoint() 00929 { 00930 BCJoint.Zero(); // basically the transformation matrix for the element 00931 BCJoint(0,1) = -1; 00932 BCJoint(0,2) = elemWidth/2; 00933 BCJoint(0,15) = 1; 00934 BCJoint(1,1) = -1; 00935 BCJoint(1,2) = -elemWidth/2; 00936 BCJoint(1,13) = 1; 00937 BCJoint(2,0) = 1; 00938 BCJoint(2,12) = -1; 00939 BCJoint(3,3) = 1; 00940 BCJoint(3,5) = elemHeight/2; 00941 BCJoint(3,12) = -1; 00942 BCJoint(4,3) = 1; 00943 BCJoint(4,5) = -elemHeight/2; 00944 BCJoint(4,14) = -1; 00945 BCJoint(5,4) = 1; 00946 BCJoint(5,13) = -1; 00947 BCJoint(6,7) = 1; 00948 BCJoint(6,8) = -elemWidth/2; 00949 BCJoint(6,15) = -1; 00950 BCJoint(7,7) = 1; 00951 BCJoint(7,8) = elemWidth/2; 00952 BCJoint(7,13) = -1; 00953 BCJoint(8,6) = 1; 00954 BCJoint(8,14) = -1; 00955 BCJoint(9,9) = -1; 00956 BCJoint(9,11) = -elemHeight/2; 00957 BCJoint(9,12) = 1; 00958 BCJoint(10,9) = -1; 00959 BCJoint(10,11) = elemHeight/2; 00960 BCJoint(10,14) = 1; 00961 BCJoint(11,10) = 1; 00962 BCJoint(11,15) = -1; 00963 BCJoint(12,12) = -1/elemHeight; 00964 BCJoint(12,13) = 1/elemWidth; 00965 BCJoint(12,14) = 1/elemHeight; 00966 BCJoint(12,15) = -1/elemWidth; 00967 00968 // based upon the new formulation 00969 BCJoint(2,2) = -(elemActHeight - elemHeight)/2; 00970 BCJoint(5,5) = -(elemActWidth - elemWidth)/2; 00971 BCJoint(8,8) = (elemActHeight - elemHeight)/2;; 00972 BCJoint(11,11) = (elemActWidth - elemWidth)/2;; 00973 00974 } 00975 00976 double BeamColumnJoint3d::getStepSize(double s0,double s1,Vector uExt,Vector duExt,Vector uInt,Vector duInt,double tol) 00977 { 00978 // finds out the factor to be used for linesearch method 00979 Vector u(16); u.Zero(); 00980 Vector fSpr(13); fSpr.Zero(); 00981 Vector kSpr(13); kSpr.Zero(); 00982 Vector intEq(4); intEq.Zero(); 00983 00984 double r0 = 0.0; // tolerance chcek for line-search 00985 double tolerance = 0.8; // slack region tolerance set for line-search 00986 00987 if (s0 != 0.0) 00988 r0 = fabs(s1/s0); 00989 00990 if (r0 <= tolerance) 00991 return 1.0; // Linsearch Not required residual decrease less than tolerance 00992 00993 if (s1 == s0) 00994 return 1.0; // Bisection would have divide by zero error if continued 00995 00996 // set some variables 00997 double etaR; 00998 double eta = 1.0; 00999 double s = s1; 01000 double etaU = 1.0; 01001 double etaL = 0.0; 01002 double sU = s1; 01003 double sL = s0; 01004 double r = r0; 01005 double etaJ = 1.0; 01006 01007 double minEta = 0.1; 01008 double maxEta = 10.0; 01009 int maxIter = 20; 01010 01011 // bracket the region 01012 int count = 0; 01013 while ((sU*sL > 0.0) && (etaU < maxEta)) { 01014 count ++; 01015 etaU *= 2.0; 01016 etaR = etaU - etaJ; 01017 01018 for (int i = 0; i < 12; i++) { 01019 u(i) = uExt(i) + duExt(i); 01020 } 01021 u(12) = uInt(0) - etaR*duInt(0); 01022 u(13) = uInt(1) - etaR*duInt(1); 01023 u(14) = uInt(2) - etaR*duInt(2); 01024 u(15) = uInt(3) - etaR*duInt(3); 01025 01026 getMatResponse(u,fSpr,kSpr); 01027 01028 intEq(0) = -fSpr(2) - fSpr(3) + fSpr(9) - fSpr(12)/elemHeight; 01029 intEq(1) = fSpr(1) - fSpr(5) - fSpr(7) + fSpr(12)/elemWidth; 01030 intEq(2) = -fSpr(4) - fSpr(8) + fSpr(10) + fSpr(12)/elemHeight; 01031 intEq(3) = fSpr(0) - fSpr(6) - fSpr(11) - fSpr(12)/elemWidth; 01032 01033 sU = duInt^intEq; 01034 01035 // check if we are happy with the solution 01036 r = fabs(sU/s0); 01037 if (r < tolerance) 01038 return etaU; 01039 01040 etaJ = etaU; 01041 } 01042 01043 if (sU*sL > 0.0) // no bracketing could be done 01044 return 1.0; 01045 01046 count = 0; 01047 while (r > tolerance && count < maxIter) { 01048 count ++; 01049 eta = (etaU + etaL)/2.0; 01050 01051 if (r > r0) eta = 1.0; 01052 01053 etaR = eta - etaJ; 01054 01055 for (int i = 0; i < 12; i++) { 01056 u(i) = uExt(i) + duExt(i); 01057 } 01058 u(12) = uInt(0) - etaR*duInt(0); 01059 u(13) = uInt(1) - etaR*duInt(1); 01060 u(14) = uInt(2) - etaR*duInt(2); 01061 u(15) = uInt(3) - etaR*duInt(3); 01062 01063 getMatResponse(u,fSpr,kSpr); 01064 01065 intEq(0) = -fSpr(2) - fSpr(3) + fSpr(9) - fSpr(12)/elemHeight; 01066 intEq(1) = fSpr(1) - fSpr(5) - fSpr(7) + fSpr(12)/elemWidth; 01067 intEq(2) = -fSpr(4) - fSpr(8) + fSpr(10) + fSpr(12)/elemHeight; 01068 intEq(3) = fSpr(0) - fSpr(6) - fSpr(11) - fSpr(12)/elemWidth; 01069 01070 s = duInt^intEq; 01071 01072 // check if we are happy with the solution 01073 r = fabs(s/s0); 01074 01075 // set variables for next iteration 01076 etaJ = eta; 01077 01078 if (s*sU < 0.0) { 01079 etaL = eta; 01080 sL = s; 01081 } else if (s*sU == 0.0) { 01082 count = maxIter; 01083 } else { 01084 etaU = eta; 01085 sU = s; 01086 } 01087 01088 if (sL == sU) 01089 count = maxIter; 01090 01091 } 01092 01093 return eta; 01094 } 01095 01096 const Matrix & 01097 BeamColumnJoint3d::getDamp(void) 01098 { 01099 // yet to do (stiffness being returned) 01100 return K; 01101 } 01102 01103 const Matrix & 01104 BeamColumnJoint3d::getMass(void) 01105 { 01106 // yet to do (stiffness being returned) 01107 return K; 01108 } 01109 01110 void 01111 BeamColumnJoint3d::zeroLoad(void) 01112 { 01113 return; 01114 } 01115 01116 int 01117 BeamColumnJoint3d::addLoad(ElementalLoad *theLoad, double loadFactor) 01118 { 01119 return 0; 01120 } 01121 01122 int 01123 BeamColumnJoint3d::addInertiaLoadToUnbalance(const Vector &accel) 01124 { 01125 return 0; 01126 } 01127 01128 01129 const Vector & 01130 BeamColumnJoint3d::getResistingForceIncInertia() 01131 { 01132 //yet to do (residual being returned) 01133 return R; 01134 } 01135 01136 int 01137 BeamColumnJoint3d::sendSelf(int commitTag, Channel &theChannel) 01138 { 01139 return -1; 01140 } 01141 01142 int 01143 BeamColumnJoint3d::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker) 01144 { 01145 return -1; 01146 } 01147 01148 int 01149 BeamColumnJoint3d::displaySelf(Renderer &theViewer, int displayMode, float fact) 01150 { 01151 return -1; // yet to do 01152 } 01153 01154 void 01155 BeamColumnJoint3d::Print(OPS_Stream &s, int flag) 01156 { 01157 s << "Element: " << this->getTag() << " Type: Beam Column Joint " << endln; 01158 for (int i = 0; i<4; i++) 01159 { 01160 s << "Node :" << connectedExternalNodes(i); 01161 s << "DOF :" << nodePtr[i]->getNumberDOF(); 01162 } 01163 s << "\nResisting Forces :" <<this->getResistingForce(); 01164 return; 01165 } 01166 01167 Response* 01168 BeamColumnJoint3d::setResponse(const char **argv, int argc, Information &eleInfo, OPS_Stream &output) 01169 { 01170 // we will compare argv[0] to determine the type of response required 01171 01172 if (strcmp(argv[0],"node1BarSlipL") == 0 || strcmp(argv[0],"node1BarslipL") ==0 || strcmp(argv[0],"Node1BarSlipL") == 0) 01173 return MaterialPtr[0]->setResponse(&argv[1], argc-1, eleInfo, output); 01174 01175 else if (strcmp(argv[0],"node1BarSlipR") == 0 || strcmp(argv[0],"node1BarslipR") ==0 || strcmp(argv[0],"Node1BarSlipR") == 0) 01176 return MaterialPtr[1]->setResponse(&argv[1], argc-1, eleInfo, output); 01177 01178 else if (strcmp(argv[0],"node1InterfaceShear") == 0 || strcmp(argv[0],"node1Interfaceshear") ==0 || strcmp(argv[0],"Node1InterfaceShear") ==0 ) 01179 return MaterialPtr[2]->setResponse(&argv[1], argc-1, eleInfo, output); 01180 01181 else if (strcmp(argv[0],"node2BarSlipB") == 0 || strcmp(argv[0],"node2BarslipB") ==0 || strcmp(argv[0],"Node2BarSlipB") == 0) 01182 return MaterialPtr[3]->setResponse(&argv[1], argc-1, eleInfo, output); 01183 01184 else if (strcmp(argv[0],"node2BarSlipT") == 0 || strcmp(argv[0],"node2BarslipT") ==0 || strcmp(argv[0],"Node2BarSlipT") == 0) 01185 return MaterialPtr[4]->setResponse(&argv[1], argc-1, eleInfo, output); 01186 01187 else if (strcmp(argv[0],"node2InterfaceShear") == 0 || strcmp(argv[0],"node2Interfaceshear") ==0 || strcmp(argv[0],"Node2InterfaceShear") ==0 ) 01188 return MaterialPtr[5]->setResponse(&argv[1], argc-1, eleInfo, output); 01189 01190 else if (strcmp(argv[0],"node3BarSlipL") == 0 || strcmp(argv[0],"node3BarslipL") ==0 || strcmp(argv[0],"Node3BarSlipL") == 0) 01191 return MaterialPtr[6]->setResponse(&argv[1], argc-1, eleInfo, output); 01192 01193 else if (strcmp(argv[0],"node3BarSlipR") == 0 || strcmp(argv[0],"node3BarslipR") ==0 || strcmp(argv[0],"Node3BarSlipR") == 0) 01194 return MaterialPtr[7]->setResponse(&argv[1], argc-1, eleInfo, output); 01195 01196 else if (strcmp(argv[0],"node3InterfaceShear") == 0 || strcmp(argv[0],"node3Interfaceshear") ==0 || strcmp(argv[0],"Node3InterfaceShear") ==0 ) 01197 return MaterialPtr[8]->setResponse(&argv[1], argc-1, eleInfo, output); 01198 01199 else if (strcmp(argv[0],"node4BarSlipB") == 0 || strcmp(argv[0],"node4BarslipB") ==0 || strcmp(argv[0],"Node4BarSlipB") == 0) 01200 return MaterialPtr[9]->setResponse(&argv[1], argc-1, eleInfo, output); 01201 01202 else if (strcmp(argv[0],"node4BarSlipT") == 0 || strcmp(argv[0],"node4BarslipT") ==0 || strcmp(argv[0],"Node4BarSlipT") == 0) 01203 return MaterialPtr[10]->setResponse(&argv[1], argc-1, eleInfo, output); 01204 01205 else if (strcmp(argv[0],"node4InterfaceShear") == 0 || strcmp(argv[0],"node4Interfaceshear") ==0 || strcmp(argv[0],"Node4InterfaceShear") ==0 ) 01206 return MaterialPtr[11]->setResponse(&argv[1], argc-1, eleInfo, output); 01207 01208 else if (strcmp(argv[0],"shearpanel") == 0 || strcmp(argv[0],"shearPanel") ==0) 01209 return MaterialPtr[12]->setResponse(&argv[1], argc-1, eleInfo, output); 01210 01211 01212 else if (strcmp(argv[0],"externalDisplacement") == 0 || strcmp(argv[0],"externaldisplacement") == 0) 01213 return new ElementResponse(this,1,Vector(12)); 01214 01215 else if (strcmp(argv[0],"internalDisplacement") == 0 || strcmp(argv[0],"internaldisplacement") == 0) 01216 return new ElementResponse(this,2,Vector(4)); 01217 01218 else if (strcmp(argv[0],"deformation") == 0 || strcmp(argv[0],"Deformation") == 0) 01219 return new ElementResponse(this,3,Vector(4)); 01220 01221 else 01222 return 0; 01223 } 01224 01225 int 01226 BeamColumnJoint3d::getResponse(int responseID, Information &eleInfo) 01227 { 01228 static Vector delta(13); 01229 static Vector def(4); 01230 static Vector U(16); 01231 static Vector Utemp(12); 01232 double bsFa, bsFb, bsFc, bsFd; 01233 double bsFac, bsFbd, isFac, isFbd; 01234 01235 switch (responseID) { 01236 case 1: 01237 if(eleInfo.theVector!=0) 01238 { 01239 (*(eleInfo.theVector))(0) = UeprCommit(0); 01240 (*(eleInfo.theVector))(1) = UeprCommit(1); 01241 (*(eleInfo.theVector))(2) = UeprCommit(2); 01242 (*(eleInfo.theVector))(3) = UeprCommit(3); 01243 (*(eleInfo.theVector))(4) = UeprCommit(4); 01244 (*(eleInfo.theVector))(5) = UeprCommit(5); 01245 (*(eleInfo.theVector))(6) = UeprCommit(6); 01246 (*(eleInfo.theVector))(7) = UeprCommit(7); 01247 (*(eleInfo.theVector))(8) = UeprCommit(8); 01248 (*(eleInfo.theVector))(9) = UeprCommit(9); 01249 (*(eleInfo.theVector))(10) = UeprCommit(10); 01250 (*(eleInfo.theVector))(11) = UeprCommit(11); 01251 (*(eleInfo.theVector))(12) = UeprCommit(12); 01252 (*(eleInfo.theVector))(13) = UeprCommit(13); 01253 (*(eleInfo.theVector))(14) = UeprCommit(14); 01254 (*(eleInfo.theVector))(15) = UeprCommit(15); 01255 (*(eleInfo.theVector))(16) = UeprCommit(16); 01256 (*(eleInfo.theVector))(17) = UeprCommit(17); 01257 (*(eleInfo.theVector))(18) = UeprCommit(18); 01258 (*(eleInfo.theVector))(19) = UeprCommit(19); 01259 (*(eleInfo.theVector))(20) = UeprCommit(20); 01260 (*(eleInfo.theVector))(21) = UeprCommit(21); 01261 (*(eleInfo.theVector))(22) = UeprCommit(22); 01262 (*(eleInfo.theVector))(23) = UeprCommit(23); 01263 } 01264 return 0; 01265 01266 case 2: 01267 if (eleInfo.theVector !=0) { 01268 (*(eleInfo.theVector))(0) = UeprIntCommit(0); 01269 (*(eleInfo.theVector))(1) = UeprIntCommit(1); 01270 (*(eleInfo.theVector))(2) = UeprIntCommit(2); 01271 (*(eleInfo.theVector))(3) = UeprIntCommit(3); 01272 } 01273 return 0; 01274 01275 case 3: 01276 01277 // modified 01.04.03 -------- determine which plane the joint lies 01278 Utemp.addMatrixVector(0.0,Transf,UeprCommit,1.0); 01279 U.Assemble(Utemp,0,1.0); 01280 U.Assemble(UeprIntCommit,12,1.0); 01281 01282 delta.addMatrixVector(0.0,BCJoint,U,1.0); 01283 01284 bsFa = fabs(delta(0) - delta(1))/elemWidth; 01285 bsFc = fabs(delta(7) - delta(6))/elemWidth; 01286 bsFac = bsFa + bsFc; 01287 bsFb = fabs(delta(4) - delta(3))/elemHeight; 01288 bsFd = fabs(delta(10) - delta(9))/elemHeight; 01289 bsFbd = bsFb + bsFd; 01290 01291 def(0) = bsFac + bsFbd; // contribution of bar slip deformation 01292 01293 01294 isFac = (delta(2) + delta(8))/elemHeight; 01295 isFbd = (delta(5) + delta(11))/elemWidth; 01296 01297 def(1) = isFac + isFbd; // contribution due to interface shear spring 01298 01299 def(2) = delta(12); // contribution due to shear panel 01300 01301 def(3) = def(0) + def(1) + def(2); // total joint deformation 01302 01303 01304 return eleInfo.setVector(def); 01305 default: 01306 return -1; 01307 } 01308 } 01309 01310 int 01311 BeamColumnJoint3d::setParameter (char **argv, int argc, Information &info) 01312 { 01313 return -1; 01314 } 01315 01316 int 01317 BeamColumnJoint3d::updateParameter (int parameterID, Information &info) 01318 { 01319 return -1; 01320 }
Generated on Mon Oct 23 15:05:08 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents