LinearCrdTransf2d.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.14 $ 00022 // $Date: 2006/09/26 18:22:24 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/coordTransformation/LinearCrdTransf2d.cpp,v $ 00024 00025 00026 // Written: Remo Magalhaes de Souza (rmsouza@ce.berkeley.edu) 00027 // Created: 04/2000 00028 // Revision: A 00029 // 00030 // Modified: May 2001 for matrix-multiply unrolling 00031 // Modified: 04/2005 Andreas Schellenberg (getBasicTrialVel, getBasicTrialAccel) 00032 // 00033 // Purpose: This file contains the implementation for the 00034 // LinearCrdTransf2d class. LinearCrdTransf2d is a linear 00035 // transformation for a planar frame between the global 00036 // and basic coordinate systems 00037 00038 00039 #include <ID.h> 00040 #include <Vector.h> 00041 #include <Matrix.h> 00042 #include <Node.h> 00043 #include <Channel.h> 00044 00045 #include <LinearCrdTransf2d.h> 00046 00047 00048 // constructor: 00049 LinearCrdTransf2d::LinearCrdTransf2d(int tag): 00050 CrdTransf2d(tag, CRDTR_TAG_LinearCrdTransf2d), 00051 nodeIPtr(0), nodeJPtr(0), 00052 nodeIOffset(0), nodeJOffset(0), 00053 cosTheta(0), sinTheta(0), L(0), 00054 nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false) 00055 { 00056 // Does nothing 00057 } 00058 00059 00060 // constructor: 00061 LinearCrdTransf2d::LinearCrdTransf2d(int tag, 00062 const Vector &rigJntOffset1, 00063 const Vector &rigJntOffset2): 00064 CrdTransf2d(tag, CRDTR_TAG_LinearCrdTransf2d), 00065 nodeIPtr(0), nodeJPtr(0), 00066 nodeIOffset(0), nodeJOffset(0), 00067 cosTheta(0), sinTheta(0), L(0), 00068 nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false) 00069 { 00070 // check rigid joint offset for node I 00071 if (&rigJntOffset1 == 0 || rigJntOffset1.Size() != 2 ) { 00072 opserr << "LinearCrdTransf2d::LinearCrdTransf2d: Invalid rigid joint offset vector for node I\n"; 00073 opserr << "Size must be 2\n"; 00074 } 00075 else if (rigJntOffset1.Norm() > 0.0) { 00076 nodeIOffset = new double[2]; 00077 nodeIOffset[0] = rigJntOffset1(0); 00078 nodeIOffset[1] = rigJntOffset1(1); 00079 } 00080 00081 // check rigid joint offset for node J 00082 if (&rigJntOffset2 == 0 || rigJntOffset2.Size() != 2 ) { 00083 opserr << "LinearCrdTransf2d::LinearCrdTransf2d: Invalid rigid joint offset vector for node J\n"; 00084 opserr << "Size must be 2\n"; 00085 } 00086 else if (rigJntOffset2.Norm() > 0.0) { 00087 nodeJOffset = new double[2]; 00088 nodeJOffset[0] = rigJntOffset2(0); 00089 nodeJOffset[1] = rigJntOffset2(1); 00090 } 00091 } 00092 00093 00094 // constructor: 00095 // invoked by a FEM_ObjectBroker, recvSelf() needs to be invoked on this object. 00096 LinearCrdTransf2d::LinearCrdTransf2d(): 00097 CrdTransf2d(0, CRDTR_TAG_LinearCrdTransf2d), 00098 nodeIPtr(0), nodeJPtr(0), 00099 nodeIOffset(0), nodeJOffset(0), 00100 cosTheta(0), sinTheta(0), L(0), 00101 nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false) 00102 { 00103 00104 } 00105 00106 00107 // destructor: 00108 LinearCrdTransf2d::~LinearCrdTransf2d() 00109 { 00110 if (nodeIOffset) 00111 delete [] nodeIOffset; 00112 if (nodeJOffset) 00113 delete [] nodeJOffset; 00114 if (nodeIInitialDisp != 0) 00115 delete [] nodeIInitialDisp; 00116 if (nodeJInitialDisp != 0) 00117 delete [] nodeJInitialDisp; 00118 } 00119 00120 00121 int 00122 LinearCrdTransf2d::commitState(void) 00123 { 00124 return 0; 00125 } 00126 00127 00128 int 00129 LinearCrdTransf2d::revertToLastCommit(void) 00130 { 00131 return 0; 00132 } 00133 00134 00135 int 00136 LinearCrdTransf2d::revertToStart(void) 00137 { 00138 return 0; 00139 } 00140 00141 00142 int 00143 LinearCrdTransf2d::initialize(Node *nodeIPointer, Node *nodeJPointer) 00144 { 00145 int error; 00146 00147 nodeIPtr = nodeIPointer; 00148 nodeJPtr = nodeJPointer; 00149 00150 if ((!nodeIPtr) || (!nodeJPtr)) 00151 { 00152 opserr << "\nLinearCrdTransf2d::initialize"; 00153 opserr << "\ninvalid pointers to the element nodes\n"; 00154 return -1; 00155 } 00156 00157 // see if there is some initial displacements at nodes 00158 if (initialDispChecked == false) { 00159 const Vector &nodeIDisp = nodeIPtr->getDisp(); 00160 const Vector &nodeJDisp = nodeJPtr->getDisp(); 00161 for (int i=0; i<3; i++) 00162 if (nodeIDisp(i) != 0.0) { 00163 nodeIInitialDisp = new double [3]; 00164 for (int j=0; j<3; j++) 00165 nodeIInitialDisp[j] = nodeIDisp(j); 00166 i = 3; 00167 } 00168 00169 for (int j=0; j<3; j++) 00170 if (nodeJDisp(j) != 0.0) { 00171 nodeJInitialDisp = new double [3]; 00172 for (int i=0; i<3; i++) 00173 nodeJInitialDisp[i] = nodeJDisp(i); 00174 j = nodeIDisp.Size(); 00175 } 00176 00177 initialDispChecked = true; 00178 } 00179 00180 // get element length and orientation 00181 if ((error = this->computeElemtLengthAndOrient())) 00182 return error; 00183 00184 return 0; 00185 } 00186 00187 00188 int 00189 LinearCrdTransf2d::update(void) 00190 { 00191 return 0; 00192 } 00193 00194 00195 int 00196 LinearCrdTransf2d::computeElemtLengthAndOrient() 00197 { 00198 // element projection 00199 static Vector dx(2); 00200 00201 const Vector &ndICoords = nodeIPtr->getCrds(); 00202 const Vector &ndJCoords = nodeJPtr->getCrds(); 00203 00204 dx(0) = ndJCoords(0) - ndICoords(0); 00205 dx(1) = ndJCoords(1) - ndICoords(1); 00206 00207 if (nodeIInitialDisp != 0) { 00208 dx(0) -= nodeIInitialDisp[0]; 00209 dx(1) -= nodeIInitialDisp[1]; 00210 } 00211 00212 if (nodeJInitialDisp != 0) { 00213 dx(0) += nodeJInitialDisp[0]; 00214 dx(1) += nodeJInitialDisp[1]; 00215 } 00216 00217 if (nodeJOffset != 0) { 00218 dx(0) += nodeJOffset[0]; 00219 dx(1) += nodeJOffset[1]; 00220 } 00221 00222 if (nodeIOffset != 0) { 00223 dx(0) -= nodeIOffset[0]; 00224 dx(1) -= nodeIOffset[1]; 00225 } 00226 00227 // calculate the element length 00228 L = dx.Norm(); 00229 00230 if (L == 0.0) 00231 { 00232 opserr << "\nLinearCrdTransf2d::computeElemtLengthAndOrien: 0 length\n"; 00233 return -2; 00234 } 00235 00236 // calculate the element local x axis components (direction cossines) 00237 // wrt to the global coordinates 00238 cosTheta = dx(0)/L; 00239 sinTheta = dx(1)/L; 00240 00241 return 0; 00242 } 00243 00244 00245 double 00246 LinearCrdTransf2d::getInitialLength(void) 00247 { 00248 return L; 00249 } 00250 00251 00252 double 00253 LinearCrdTransf2d::getDeformedLength(void) 00254 { 00255 return L; 00256 } 00257 00258 00259 const Vector & 00260 LinearCrdTransf2d::getBasicTrialDisp (void) 00261 { 00262 // determine global displacements 00263 const Vector &disp1 = nodeIPtr->getTrialDisp(); 00264 const Vector &disp2 = nodeJPtr->getTrialDisp(); 00265 00266 static double ug[6]; 00267 for (int i = 0; i < 3; i++) { 00268 ug[i] = disp1(i); 00269 ug[i+3] = disp2(i); 00270 } 00271 00272 if (nodeIInitialDisp != 0) { 00273 for (int j=0; j<3; j++) 00274 ug[j] -= nodeIInitialDisp[j]; 00275 } 00276 00277 if (nodeJInitialDisp != 0) { 00278 for (int j=0; j<3; j++) 00279 ug[j+3] -= nodeJInitialDisp[j]; 00280 } 00281 00282 static Vector ub(3); 00283 00284 double oneOverL = 1.0/L; 00285 double sl = sinTheta*oneOverL; 00286 double cl = cosTheta*oneOverL; 00287 00288 ub(0) = -cosTheta*ug[0] - sinTheta*ug[1] + 00289 cosTheta*ug[3] + sinTheta*ug[4]; 00290 00291 ub(1) = -sl*ug[0] + cl*ug[1] + ug[2] + 00292 sl*ug[3] - cl*ug[4]; 00293 00294 if (nodeIOffset != 0) { 00295 double t02 = -cosTheta*nodeIOffset[1] + sinTheta*nodeIOffset[0]; 00296 double t12 = sinTheta*nodeIOffset[1] + cosTheta*nodeIOffset[0]; 00297 ub(0) -= t02*ug[2]; 00298 ub(1) += oneOverL*t12*ug[2]; 00299 } 00300 00301 if (nodeJOffset != 0) { 00302 double t35 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00303 double t45 = sinTheta*nodeJOffset[1] + cosTheta*nodeJOffset[0]; 00304 ub(0) += t35*ug[5]; 00305 ub(1) -= oneOverL*t45*ug[5]; 00306 } 00307 00308 ub(2) = ub(1) + ug[5] - ug[2]; 00309 00310 return ub; 00311 } 00312 00313 00314 const Vector & 00315 LinearCrdTransf2d::getBasicIncrDisp (void) 00316 { 00317 // determine global displacements 00318 const Vector &disp1 = nodeIPtr->getIncrDisp(); 00319 const Vector &disp2 = nodeJPtr->getIncrDisp(); 00320 00321 static double dug[6]; 00322 for (int i = 0; i < 3; i++) { 00323 dug[i] = disp1(i); 00324 dug[i+3] = disp2(i); 00325 } 00326 00327 static Vector dub(3); 00328 00329 double oneOverL = 1.0/L; 00330 double sl = sinTheta*oneOverL; 00331 double cl = cosTheta*oneOverL; 00332 00333 dub(0) = -cosTheta*dug[0] - sinTheta*dug[1] + 00334 cosTheta*dug[3] + sinTheta*dug[4]; 00335 00336 dub(1) = -sl*dug[0] + cl*dug[1] + dug[2] + 00337 sl*dug[3] - cl*dug[4]; 00338 00339 if (nodeIOffset != 0) { 00340 double t02 = -cosTheta*nodeIOffset[1] + sinTheta*nodeIOffset[0]; 00341 double t12 = sinTheta*nodeIOffset[1] + cosTheta*nodeIOffset[0]; 00342 dub(0) -= t02*dug[2]; 00343 dub(1) += oneOverL*t12*dug[2]; 00344 } 00345 00346 if (nodeJOffset != 0) { 00347 double t35 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00348 double t45 = sinTheta*nodeJOffset[1] + cosTheta*nodeJOffset[0]; 00349 dub(0) += t35*dug[5]; 00350 dub(1) -= oneOverL*t45*dug[5]; 00351 } 00352 00353 dub(2) = dub(1) + dug[5] - dug[2]; 00354 00355 return dub; 00356 } 00357 00358 00359 const Vector & 00360 LinearCrdTransf2d::getBasicIncrDeltaDisp(void) 00361 { 00362 // determine global displacements 00363 const Vector &disp1 = nodeIPtr->getIncrDeltaDisp(); 00364 const Vector &disp2 = nodeJPtr->getIncrDeltaDisp(); 00365 00366 static double Dug[6]; 00367 for (int i = 0; i < 3; i++) { 00368 Dug[i] = disp1(i); 00369 Dug[i+3] = disp2(i); 00370 } 00371 00372 static Vector Dub(3); 00373 00374 double oneOverL = 1.0/L; 00375 double sl = sinTheta*oneOverL; 00376 double cl = cosTheta*oneOverL; 00377 00378 Dub(0) = -cosTheta*Dug[0] - sinTheta*Dug[1] + 00379 cosTheta*Dug[3] + sinTheta*Dug[4]; 00380 00381 Dub(1) = -sl*Dug[0] + cl*Dug[1] + Dug[2] + 00382 sl*Dug[3] - cl*Dug[4]; 00383 00384 if (nodeIOffset != 0) { 00385 double t02 = -cosTheta*nodeIOffset[1] + sinTheta*nodeIOffset[0]; 00386 double t12 = sinTheta*nodeIOffset[1] + cosTheta*nodeIOffset[0]; 00387 Dub(0) -= t02*Dug[2]; 00388 Dub(1) += oneOverL*t12*Dug[2]; 00389 } 00390 00391 if (nodeJOffset != 0) { 00392 double t35 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00393 double t45 = sinTheta*nodeJOffset[1] + cosTheta*nodeJOffset[0]; 00394 Dub(0) += t35*Dug[5]; 00395 Dub(1) -= oneOverL*t45*Dug[5]; 00396 } 00397 00398 Dub(2) = Dub(1) + Dug[5] - Dug[2]; 00399 00400 return Dub; 00401 } 00402 00403 00404 const Vector & 00405 LinearCrdTransf2d::getBasicTrialVel(void) 00406 { 00407 // determine global velocities 00408 const Vector &vel1 = nodeIPtr->getTrialVel(); 00409 const Vector &vel2 = nodeJPtr->getTrialVel(); 00410 00411 static double vg[6]; 00412 for (int i = 0; i < 3; i++) { 00413 vg[i] = vel1(i); 00414 vg[i+3] = vel2(i); 00415 } 00416 00417 static Vector vb(3); 00418 00419 double oneOverL = 1.0/L; 00420 double sl = sinTheta*oneOverL; 00421 double cl = cosTheta*oneOverL; 00422 00423 vb(0) = -cosTheta*vg[0] - sinTheta*vg[1] + 00424 cosTheta*vg[3] + sinTheta*vg[4]; 00425 00426 vb(1) = -sl*vg[0] + cl*vg[1] + vg[2] + 00427 sl*vg[3] - cl*vg[4]; 00428 00429 if (nodeIOffset != 0) { 00430 double t02 = -cosTheta*nodeIOffset[1] + sinTheta*nodeIOffset[0]; 00431 double t12 = sinTheta*nodeIOffset[1] + cosTheta*nodeIOffset[0]; 00432 vb(0) -= t02*vg[2]; 00433 vb(1) += oneOverL*t12*vg[2]; 00434 } 00435 00436 if (nodeJOffset != 0) { 00437 double t35 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00438 double t45 = sinTheta*nodeJOffset[1] + cosTheta*nodeJOffset[0]; 00439 vb(0) += t35*vg[5]; 00440 vb(1) -= oneOverL*t45*vg[5]; 00441 } 00442 00443 vb(2) = vb(1) + vg[5] - vg[2]; 00444 00445 return vb; 00446 } 00447 00448 00449 const Vector & 00450 LinearCrdTransf2d::getBasicTrialAccel(void) 00451 { 00452 // determine global accelerations 00453 const Vector &accel1 = nodeIPtr->getTrialAccel(); 00454 const Vector &accel2 = nodeJPtr->getTrialAccel(); 00455 00456 static double ag[6]; 00457 for (int i = 0; i < 3; i++) { 00458 ag[i] = accel1(i); 00459 ag[i+3] = accel2(i); 00460 } 00461 00462 static Vector ab(3); 00463 00464 double oneOverL = 1.0/L; 00465 double sl = sinTheta*oneOverL; 00466 double cl = cosTheta*oneOverL; 00467 00468 ab(0) = -cosTheta*ag[0] - sinTheta*ag[1] + 00469 cosTheta*ag[3] + sinTheta*ag[4]; 00470 00471 ab(1) = -sl*ag[0] + cl*ag[1] + ag[2] + 00472 sl*ag[3] - cl*ag[4]; 00473 00474 if (nodeIOffset != 0) { 00475 double t02 = -cosTheta*nodeIOffset[1] + sinTheta*nodeIOffset[0]; 00476 double t12 = sinTheta*nodeIOffset[1] + cosTheta*nodeIOffset[0]; 00477 ab(0) -= t02*ag[2]; 00478 ab(1) += oneOverL*t12*ag[2]; 00479 } 00480 00481 if (nodeJOffset != 0) { 00482 double t35 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00483 double t45 = sinTheta*nodeJOffset[1] + cosTheta*nodeJOffset[0]; 00484 ab(0) += t35*ag[5]; 00485 ab(1) -= oneOverL*t45*ag[5]; 00486 } 00487 00488 ab(2) = ab(1) + ag[5] - ag[2]; 00489 00490 return ab; 00491 } 00492 00493 00494 const Vector & 00495 LinearCrdTransf2d::getBasicTrialDispShapeSensitivity (void) 00496 { 00497 // Want to return dAdh * u 00498 00499 // determine global displacements 00500 const Vector &disp1 = nodeIPtr->getTrialDisp(); 00501 const Vector &disp2 = nodeJPtr->getTrialDisp(); 00502 00503 static double ug[6]; 00504 for (int i = 0; i < 3; i++) { 00505 ug[i] = disp1(i); 00506 ug[i+3] = disp2(i); 00507 } 00508 00509 if (nodeIInitialDisp != 0) { 00510 for (int j=0; j<3; j++) 00511 ug[j] -= nodeIInitialDisp[j]; 00512 } 00513 00514 if (nodeJInitialDisp != 0) { 00515 for (int j=0; j<3; j++) 00516 ug[j+3] -= nodeJInitialDisp[j]; 00517 } 00518 00519 static Vector ub(3); 00520 ub.Zero(); 00521 00522 static ID nodeParameterID(2); 00523 nodeParameterID(0) = nodeIPtr->getCrdsSensitivity(); 00524 nodeParameterID(1) = nodeJPtr->getCrdsSensitivity(); 00525 00526 if (nodeParameterID(0) != 0 || nodeParameterID(1) != 0) { 00527 00528 if (nodeIOffset != 0 || nodeJOffset != 0) { 00529 opserr << "ERROR: Currently a node offset cannot be used in " << endln 00530 << " conjunction with random nodal coordinates." << endln; 00531 } 00532 00533 double dcosdh, dsindh, dsldh, dcldh; 00534 00535 double dx = cosTheta*L; 00536 double dy = sinTheta*L; 00537 00538 if (nodeParameterID(0) == 1) { // here x1 is random 00539 dcosdh = (-L+dx*dx/L)/(L*L); 00540 dsindh = dx*dy/(L*L*L); 00541 dcldh = (-L*L+dx*dx*2)/(L*L*L*L); 00542 dsldh = 2*dx*dy/(L*L*L*L); 00543 } 00544 if (nodeParameterID(0) == 2) { // here y1 is random 00545 dsindh = (-L+dy*dy/L)/(L*L); 00546 dcosdh = dx*dy/(L*L*L); 00547 dsldh = (-L*L+dy*dy*2)/(L*L*L*L); 00548 dcldh = 2*dx*dy/(L*L*L*L); 00549 } 00550 00551 if (nodeParameterID(1) == 1) { // here x2 is random 00552 dcosdh = (L-dx*dx/L)/(L*L); 00553 dsindh = -dx*dy/(L*L*L); 00554 dcldh = (L*L-dx*dx*2)/(L*L*L*L); 00555 dsldh = -2*dx*dy/(L*L*L*L); 00556 } 00557 if (nodeParameterID(1) == 2) { // here y2 is random 00558 dsindh = (L-dy*dy/L)/(L*L); 00559 dcosdh = -dx*dy/(L*L*L); 00560 dsldh = (L*L-dy*dy*2)/(L*L*L*L); 00561 dcldh = -2*dx*dy/(L*L*L*L); 00562 } 00563 00564 ub(0) = -dcosdh*ug[0] - dsindh*ug[1] + dcosdh*ug[3] + dsindh*ug[4]; 00565 00566 ub(1) = -dsldh*ug[0] + dcldh*ug[1] + dsldh*ug[3] - dcldh*ug[4]; 00567 00568 ub(2) = ub(1); 00569 } 00570 00571 return ub; 00572 } 00573 00574 00575 const Vector & 00576 LinearCrdTransf2d::getGlobalResistingForce(const Vector &pb, const Vector &p0) 00577 { 00578 // transform resisting forces from the basic system to local coordinates 00579 static double pl[6]; 00580 00581 double q0 = pb(0); 00582 double q1 = pb(1); 00583 double q2 = pb(2); 00584 00585 double oneOverL = 1.0/L; 00586 00587 double V = oneOverL*(q1+q2); 00588 pl[0] = -q0; 00589 pl[1] = V; 00590 pl[2] = q1; 00591 pl[3] = q0; 00592 pl[4] = -V; 00593 pl[5] = q2; 00594 00595 // add end forces due to element p0 loads 00596 pl[0] += p0(0); 00597 pl[1] += p0(1); 00598 pl[4] += p0(2); 00599 00600 // transform resisting forces from local to global coordinates 00601 static Vector pg(6); 00602 00603 pg(0) = cosTheta*pl[0] - sinTheta*pl[1]; 00604 pg(1) = sinTheta*pl[0] + cosTheta*pl[1]; 00605 00606 pg(3) = cosTheta*pl[3] - sinTheta*pl[4]; 00607 pg(4) = sinTheta*pl[3] + cosTheta*pl[4]; 00608 00609 pg(2) = pl[2]; 00610 pg(5) = pl[5]; 00611 00612 if (nodeIOffset != 0) { 00613 double t02 = -cosTheta*nodeIOffset[1] + sinTheta*nodeIOffset[0]; 00614 double t12 = sinTheta*nodeIOffset[1] + cosTheta*nodeIOffset[0]; 00615 00616 pg(2) += t02*pl[0] + t12*pl[1]; 00617 } 00618 00619 if (nodeJOffset != 0) { 00620 double t35 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00621 double t45 = sinTheta*nodeJOffset[1] + cosTheta*nodeJOffset[0]; 00622 00623 pg(5) += t35*pl[3] + t45*pl[4]; 00624 } 00625 00626 return pg; 00627 } 00628 00629 00630 const Matrix & 00631 LinearCrdTransf2d::getGlobalStiffMatrix (const Matrix &kb, const Vector &pb) 00632 { 00633 static Matrix kg(6,6); 00634 static double tmp [6][6]; 00635 00636 double oneOverL = 1.0/L; 00637 00638 double kb00, kb01, kb02, kb10, kb11, kb12, kb20, kb21, kb22; 00639 00640 kb00 = kb(0,0); kb01 = kb(0,1); kb02 = kb(0,2); 00641 kb10 = kb(1,0); kb11 = kb(1,1); kb12 = kb(1,2); 00642 kb20 = kb(2,0); kb21 = kb(2,1); kb22 = kb(2,2); 00643 00644 double t02 = 0.0; 00645 double t12 = 1.0; 00646 double t22 = 0.0; 00647 00648 if (nodeIOffset != 0) { 00649 t02 = cosTheta*nodeIOffset[1] - sinTheta*nodeIOffset[0]; 00650 t12 = oneOverL*(sinTheta*nodeIOffset[1]+cosTheta*nodeIOffset[0]) + 1.0; 00651 t22 = oneOverL*(sinTheta*nodeIOffset[1]+cosTheta*nodeIOffset[0]); 00652 } 00653 00654 double t05 = 0.0; 00655 double t15 = 0.0; 00656 double t25 = 1.0; 00657 00658 if (nodeJOffset != 0) { 00659 t05 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00660 t15 = -oneOverL*(sinTheta*nodeJOffset[1]+cosTheta*nodeJOffset[0]); 00661 t25 = -oneOverL*(sinTheta*nodeJOffset[1]+cosTheta*nodeJOffset[0]) + 1.0; 00662 } 00663 00664 double sl = sinTheta*oneOverL; 00665 double cl = cosTheta*oneOverL; 00666 00667 tmp[0][0] = -cosTheta*kb00 - sl*(kb01+kb02); 00668 tmp[0][1] = -sinTheta*kb00 + cl*(kb01+kb02); 00669 tmp[0][2] = (nodeIOffset) ? t02*kb00 + t12*kb01 + t22*kb02 : kb01; 00670 tmp[0][3] = -tmp[0][0]; 00671 tmp[0][4] = -tmp[0][1]; 00672 tmp[0][5] = (nodeJOffset) ? t05*kb00 + t15*kb01 + t25*kb02 : kb02; 00673 00674 tmp[1][0] = -cosTheta*kb10 - sl*(kb11+kb12); 00675 tmp[1][1] = -sinTheta*kb10 + cl*(kb11+kb12); 00676 tmp[1][2] = (nodeIOffset) ? t02*kb10 + t12*kb11 + t22*kb12 : kb11; 00677 tmp[1][3] = -tmp[1][0]; 00678 tmp[1][4] = -tmp[1][1]; 00679 tmp[1][5] = (nodeJOffset) ? t05*kb10 + t15*kb11 + t25*kb12 : kb12; 00680 00681 tmp[2][0] = -cosTheta*kb20 - sl*(kb21+kb22); 00682 tmp[2][1] = -sinTheta*kb20 + cl*(kb21+kb22); 00683 tmp[2][2] = (nodeIOffset) ? t02*kb20 + t12*kb21 + t22*kb22 : kb21; 00684 tmp[2][3] = -tmp[2][0]; 00685 tmp[2][4] = -tmp[2][1]; 00686 tmp[2][5] = (nodeJOffset) ? t05*kb20 + t15*kb21 + t25*kb22 : kb22; 00687 00688 kg(0,0) = -cosTheta*tmp[0][0] - sl*(tmp[1][0]+tmp[2][0]); 00689 kg(0,1) = -cosTheta*tmp[0][1] - sl*(tmp[1][1]+tmp[2][1]); 00690 kg(0,2) = -cosTheta*tmp[0][2] - sl*(tmp[1][2]+tmp[2][2]); 00691 kg(0,3) = -cosTheta*tmp[0][3] - sl*(tmp[1][3]+tmp[2][3]); 00692 kg(0,4) = -cosTheta*tmp[0][4] - sl*(tmp[1][4]+tmp[2][4]); 00693 kg(0,5) = -cosTheta*tmp[0][5] - sl*(tmp[1][5]+tmp[2][5]); 00694 00695 kg(1,0) = -sinTheta*tmp[0][0] + cl*(tmp[1][0]+tmp[2][0]); 00696 kg(1,1) = -sinTheta*tmp[0][1] + cl*(tmp[1][1]+tmp[2][1]); 00697 kg(1,2) = -sinTheta*tmp[0][2] + cl*(tmp[1][2]+tmp[2][2]); 00698 kg(1,3) = -sinTheta*tmp[0][3] + cl*(tmp[1][3]+tmp[2][3]); 00699 kg(1,4) = -sinTheta*tmp[0][4] + cl*(tmp[1][4]+tmp[2][4]); 00700 kg(1,5) = -sinTheta*tmp[0][5] + cl*(tmp[1][5]+tmp[2][5]); 00701 00702 if (nodeIOffset) { 00703 kg(2,0) = t02*tmp[0][0] + t12*tmp[1][0] + t22*tmp[2][0]; 00704 kg(2,1) = t02*tmp[0][1] + t12*tmp[1][1] + t22*tmp[2][1]; 00705 kg(2,2) = t02*tmp[0][2] + t12*tmp[1][2] + t22*tmp[2][2]; 00706 kg(2,3) = t02*tmp[0][3] + t12*tmp[1][3] + t22*tmp[2][3]; 00707 kg(2,4) = t02*tmp[0][4] + t12*tmp[1][4] + t22*tmp[2][4]; 00708 kg(2,5) = t02*tmp[0][5] + t12*tmp[1][5] + t22*tmp[2][5]; 00709 } 00710 else { 00711 kg(2,0) = tmp[1][0]; 00712 kg(2,1) = tmp[1][1]; 00713 kg(2,2) = tmp[1][2]; 00714 kg(2,3) = tmp[1][3]; 00715 kg(2,4) = tmp[1][4]; 00716 kg(2,5) = tmp[1][5]; 00717 } 00718 00719 kg(3,0) = -kg(0,0); 00720 kg(3,1) = -kg(0,1); 00721 kg(3,2) = -kg(0,2); 00722 kg(3,3) = -kg(0,3); 00723 kg(3,4) = -kg(0,4); 00724 kg(3,5) = -kg(0,5); 00725 00726 kg(4,0) = -kg(1,0); 00727 kg(4,1) = -kg(1,1); 00728 kg(4,2) = -kg(1,2); 00729 kg(4,3) = -kg(1,3); 00730 kg(4,4) = -kg(1,4); 00731 kg(4,5) = -kg(1,5); 00732 00733 if (nodeJOffset) { 00734 kg(5,0) = t05*tmp[0][0] + t15*tmp[1][0] + t25*tmp[2][0]; 00735 kg(5,1) = t05*tmp[0][1] + t15*tmp[1][1] + t25*tmp[2][1]; 00736 kg(5,2) = t05*tmp[0][2] + t15*tmp[1][2] + t25*tmp[2][2]; 00737 kg(5,3) = t05*tmp[0][3] + t15*tmp[1][3] + t25*tmp[2][3]; 00738 kg(5,4) = t05*tmp[0][4] + t15*tmp[1][4] + t25*tmp[2][4]; 00739 kg(5,5) = t05*tmp[0][5] + t15*tmp[1][5] + t25*tmp[2][5]; 00740 } 00741 else { 00742 kg(5,0) = tmp[2][0]; 00743 kg(5,1) = tmp[2][1]; 00744 kg(5,2) = tmp[2][2]; 00745 kg(5,3) = tmp[2][3]; 00746 kg(5,4) = tmp[2][4]; 00747 kg(5,5) = tmp[2][5]; 00748 } 00749 00750 return kg; 00751 } 00752 00753 00754 const Matrix & 00755 LinearCrdTransf2d::getInitialGlobalStiffMatrix (const Matrix &kb) 00756 { 00757 static Matrix kg(6,6); 00758 static double tmp [6][6]; 00759 00760 double oneOverL = 1.0/L; 00761 00762 double kb00, kb01, kb02, kb10, kb11, kb12, kb20, kb21, kb22; 00763 00764 kb00 = kb(0,0); kb01 = kb(0,1); kb02 = kb(0,2); 00765 kb10 = kb(1,0); kb11 = kb(1,1); kb12 = kb(1,2); 00766 kb20 = kb(2,0); kb21 = kb(2,1); kb22 = kb(2,2); 00767 00768 double t02 = 0.0; 00769 double t12 = 1.0; 00770 double t22 = 0.0; 00771 00772 if (nodeIOffset != 0) { 00773 t02 = cosTheta*nodeIOffset[1] - sinTheta*nodeIOffset[0]; 00774 t12 = oneOverL*(sinTheta*nodeIOffset[1]+cosTheta*nodeIOffset[0]) + 1.0; 00775 t22 = oneOverL*(sinTheta*nodeIOffset[1]+cosTheta*nodeIOffset[0]); 00776 } 00777 00778 double t05 = 0.0; 00779 double t15 = 0.0; 00780 double t25 = 1.0; 00781 00782 if (nodeJOffset != 0) { 00783 t05 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 00784 t15 = -oneOverL*(sinTheta*nodeJOffset[1]+cosTheta*nodeJOffset[0]); 00785 t25 = -oneOverL*(sinTheta*nodeJOffset[1]+cosTheta*nodeJOffset[0]) + 1.0; 00786 } 00787 00788 double sl = sinTheta*oneOverL; 00789 double cl = cosTheta*oneOverL; 00790 00791 tmp[0][0] = -cosTheta*kb00 - sl*(kb01+kb02); 00792 tmp[0][1] = -sinTheta*kb00 + cl*(kb01+kb02); 00793 tmp[0][2] = (nodeIOffset) ? t02*kb00 + t12*kb01 + t22*kb02 : kb01; 00794 tmp[0][3] = -tmp[0][0]; 00795 tmp[0][4] = -tmp[0][1]; 00796 tmp[0][5] = (nodeJOffset) ? t05*kb00 + t15*kb01 + t25*kb02 : kb02; 00797 00798 tmp[1][0] = -cosTheta*kb10 - sl*(kb11+kb12); 00799 tmp[1][1] = -sinTheta*kb10 + cl*(kb11+kb12); 00800 tmp[1][2] = (nodeIOffset) ? t02*kb10 + t12*kb11 + t22*kb12 : kb11; 00801 tmp[1][3] = -tmp[1][0]; 00802 tmp[1][4] = -tmp[1][1]; 00803 tmp[1][5] = (nodeJOffset) ? t05*kb10 + t15*kb11 + t25*kb12 : kb12; 00804 00805 tmp[2][0] = -cosTheta*kb20 - sl*(kb21+kb22); 00806 tmp[2][1] = -sinTheta*kb20 + cl*(kb21+kb22); 00807 tmp[2][2] = (nodeIOffset) ? t02*kb20 + t12*kb21 + t22*kb22 : kb21; 00808 tmp[2][3] = -tmp[2][0]; 00809 tmp[2][4] = -tmp[2][1]; 00810 tmp[2][5] = (nodeJOffset) ? t05*kb20 + t15*kb21 + t25*kb22 : kb22; 00811 00812 kg(0,0) = -cosTheta*tmp[0][0] - sl*(tmp[1][0]+tmp[2][0]); 00813 kg(0,1) = -cosTheta*tmp[0][1] - sl*(tmp[1][1]+tmp[2][1]); 00814 kg(0,2) = -cosTheta*tmp[0][2] - sl*(tmp[1][2]+tmp[2][2]); 00815 kg(0,3) = -cosTheta*tmp[0][3] - sl*(tmp[1][3]+tmp[2][3]); 00816 kg(0,4) = -cosTheta*tmp[0][4] - sl*(tmp[1][4]+tmp[2][4]); 00817 kg(0,5) = -cosTheta*tmp[0][5] - sl*(tmp[1][5]+tmp[2][5]); 00818 00819 kg(1,0) = -sinTheta*tmp[0][0] + cl*(tmp[1][0]+tmp[2][0]); 00820 kg(1,1) = -sinTheta*tmp[0][1] + cl*(tmp[1][1]+tmp[2][1]); 00821 kg(1,2) = -sinTheta*tmp[0][2] + cl*(tmp[1][2]+tmp[2][2]); 00822 kg(1,3) = -sinTheta*tmp[0][3] + cl*(tmp[1][3]+tmp[2][3]); 00823 kg(1,4) = -sinTheta*tmp[0][4] + cl*(tmp[1][4]+tmp[2][4]); 00824 kg(1,5) = -sinTheta*tmp[0][5] + cl*(tmp[1][5]+tmp[2][5]); 00825 00826 if (nodeIOffset) { 00827 kg(2,0) = t02*tmp[0][0] + t12*tmp[1][0] + t22*tmp[2][0]; 00828 kg(2,1) = t02*tmp[0][1] + t12*tmp[1][1] + t22*tmp[2][1]; 00829 kg(2,2) = t02*tmp[0][2] + t12*tmp[1][2] + t22*tmp[2][2]; 00830 kg(2,3) = t02*tmp[0][3] + t12*tmp[1][3] + t22*tmp[2][3]; 00831 kg(2,4) = t02*tmp[0][4] + t12*tmp[1][4] + t22*tmp[2][4]; 00832 kg(2,5) = t02*tmp[0][5] + t12*tmp[1][5] + t22*tmp[2][5]; 00833 } 00834 else { 00835 kg(2,0) = tmp[1][0]; 00836 kg(2,1) = tmp[1][1]; 00837 kg(2,2) = tmp[1][2]; 00838 kg(2,3) = tmp[1][3]; 00839 kg(2,4) = tmp[1][4]; 00840 kg(2,5) = tmp[1][5]; 00841 } 00842 00843 kg(3,0) = -kg(0,0); 00844 kg(3,1) = -kg(0,1); 00845 kg(3,2) = -kg(0,2); 00846 kg(3,3) = -kg(0,3); 00847 kg(3,4) = -kg(0,4); 00848 kg(3,5) = -kg(0,5); 00849 00850 kg(4,0) = -kg(1,0); 00851 kg(4,1) = -kg(1,1); 00852 kg(4,2) = -kg(1,2); 00853 kg(4,3) = -kg(1,3); 00854 kg(4,4) = -kg(1,4); 00855 kg(4,5) = -kg(1,5); 00856 00857 if (nodeJOffset) { 00858 kg(5,0) = t05*tmp[0][0] + t15*tmp[1][0] + t25*tmp[2][0]; 00859 kg(5,1) = t05*tmp[0][1] + t15*tmp[1][1] + t25*tmp[2][1]; 00860 kg(5,2) = t05*tmp[0][2] + t15*tmp[1][2] + t25*tmp[2][2]; 00861 kg(5,3) = t05*tmp[0][3] + t15*tmp[1][3] + t25*tmp[2][3]; 00862 kg(5,4) = t05*tmp[0][4] + t15*tmp[1][4] + t25*tmp[2][4]; 00863 kg(5,5) = t05*tmp[0][5] + t15*tmp[1][5] + t25*tmp[2][5]; 00864 } 00865 else { 00866 kg(5,0) = tmp[2][0]; 00867 kg(5,1) = tmp[2][1]; 00868 kg(5,2) = tmp[2][2]; 00869 kg(5,3) = tmp[2][3]; 00870 kg(5,4) = tmp[2][4]; 00871 kg(5,5) = tmp[2][5]; 00872 } 00873 00874 return kg; 00875 } 00876 00877 00878 CrdTransf2d * 00879 LinearCrdTransf2d::getCopy(void) 00880 { 00881 // create a new instance of LinearCrdTransf2d 00882 00883 LinearCrdTransf2d *theCopy; 00884 00885 Vector offsetI(2); 00886 Vector offsetJ(2); 00887 00888 if (nodeIOffset != 0) { 00889 offsetI(0) = nodeIOffset[0]; 00890 offsetI(1) = nodeIOffset[1]; 00891 } 00892 00893 if (nodeJOffset != 0) { 00894 offsetJ(0) = nodeJOffset[0]; 00895 offsetJ(1) = nodeJOffset[1]; 00896 } 00897 00898 theCopy = new LinearCrdTransf2d(this->getTag(), offsetI, offsetJ); 00899 00900 theCopy->nodeIPtr = nodeIPtr; 00901 theCopy->nodeJPtr = nodeJPtr; 00902 theCopy->cosTheta = cosTheta; 00903 theCopy->sinTheta = sinTheta; 00904 theCopy->L = L; 00905 00906 return theCopy; 00907 } 00908 00909 00910 int 00911 LinearCrdTransf2d::sendSelf(int cTag, Channel &theChannel) 00912 { 00913 int res = 0; 00914 00915 static Vector data(12); 00916 data(0) = this->getTag(); 00917 data(1) = L; 00918 if (nodeIOffset != 0) { 00919 data(2) = nodeIOffset[0]; 00920 data(3) = nodeIOffset[1]; 00921 } else { 00922 data(2) = 0.0; 00923 data(3) = 0.0; 00924 } 00925 00926 if (nodeJOffset != 0) { 00927 data(4) = nodeJOffset[0]; 00928 data(5) = nodeJOffset[1]; 00929 } else { 00930 data(4) = 0.0; 00931 data(5) = 0.0; 00932 } 00933 00934 if (nodeIInitialDisp != 0) { 00935 data(6) = nodeIInitialDisp[0]; 00936 data(7) = nodeIInitialDisp[1]; 00937 data(8) = nodeIInitialDisp[2]; 00938 } else { 00939 data(6) = 0.0; 00940 data(7) = 0.0; 00941 data(8) = 0.0; 00942 } 00943 00944 if (nodeJInitialDisp != 0) { 00945 data(9) = nodeJInitialDisp[0]; 00946 data(10) = nodeJInitialDisp[1]; 00947 data(11) = nodeJInitialDisp[2]; 00948 } else { 00949 data(9) = 0.0; 00950 data(10) = 0.0; 00951 data(11) = 0.0; 00952 } 00953 00954 res += theChannel.sendVector(this->getDbTag(), cTag, data); 00955 if (res < 0) { 00956 opserr << "LinearCrdTransf2d::sendSelf - failed to send Vector\n"; 00957 return res; 00958 } 00959 00960 return res; 00961 } 00962 00963 00964 int 00965 LinearCrdTransf2d::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker) 00966 { 00967 int res = 0; 00968 00969 static Vector data(12); 00970 00971 res += theChannel.recvVector(this->getDbTag(), cTag, data); 00972 if (res < 0) { 00973 opserr << "LinearCrdTransf2d::recvSelf - failed to receive Vector\n"; 00974 return res; 00975 } 00976 00977 this->setTag((int)data(0)); 00978 L = data(1); 00979 data(0) = this->getTag(); 00980 data(1) = L; 00981 00982 int flag; 00983 int i,j; 00984 00985 flag = 0; 00986 for (i=2; i<=3; i++) 00987 if (data(i) != 0.0) 00988 flag = 1; 00989 if (flag == 1) { 00990 if (nodeIOffset == 0) 00991 nodeIOffset = new double[2]; 00992 for (i=2, j=0; i<=3; i++, j++) 00993 nodeIOffset[j] = data(i); 00994 } 00995 00996 flag = 0; 00997 for (i=4; i<=5; i++) 00998 if (data(i) != 0.0) 00999 flag = 1; 01000 if (flag == 1) { 01001 if (nodeJOffset == 0) 01002 nodeJOffset = new double[2]; 01003 for (i=4, j=0; i<=5; i++, j++) 01004 nodeJOffset[j] = data(i); 01005 } 01006 01007 flag = 0; 01008 for (i=6; i<=8; i++) 01009 if (data(i) != 0.0) 01010 flag = 1; 01011 if (flag == 1) { 01012 if (nodeIInitialDisp == 0) 01013 nodeIInitialDisp = new double[3]; 01014 for (i=6, j=0; i<=7; i++, j++) 01015 nodeIInitialDisp[j] = data(i); 01016 } 01017 01018 flag = 0; 01019 for (i=9; i<=11; i++) 01020 if (data(i) != 0.0) 01021 flag = 1; 01022 if (flag == 1) { 01023 if (nodeJInitialDisp == 0) 01024 nodeJInitialDisp = new double [3]; 01025 for (i=9, j=0; i<=11; i++, j++) 01026 nodeJInitialDisp[j] = data(i); 01027 } 01028 01029 initialDispChecked = true; 01030 01031 return res; 01032 } 01033 01034 01035 const Vector & 01036 LinearCrdTransf2d::getPointGlobalCoordFromLocal(const Vector &xl) 01037 { 01038 static Vector xg(2); 01039 01040 const Vector &nodeICoords = nodeIPtr->getCrds(); 01041 xg(0) = nodeICoords(0); 01042 xg(1) = nodeICoords(1); 01043 01044 if (nodeIOffset) { 01045 xg(0) += nodeIOffset[0]; 01046 xg(1) += nodeIOffset[1]; 01047 } 01048 01049 // xg = xg + Rlj'*xl 01050 xg(0) += cosTheta*xl(0) - sinTheta*xl(1); 01051 xg(1) += sinTheta*xl(0) + cosTheta*xl(1); 01052 01053 return xg; 01054 } 01055 01056 01057 const Vector & 01058 LinearCrdTransf2d::getPointGlobalDisplFromBasic (double xi, const Vector &uxb) 01059 { 01060 // determine global displacements 01061 const Vector &disp1 = nodeIPtr->getTrialDisp(); 01062 const Vector &disp2 = nodeJPtr->getTrialDisp(); 01063 01064 static Vector ug(6); 01065 for (int i = 0; i < 3; i++) 01066 { 01067 ug(i) = disp1(i); 01068 ug(i+3) = disp2(i); 01069 } 01070 01071 if (nodeIInitialDisp != 0) { 01072 for (int j=0; j<3; j++) 01073 ug[j] -= nodeIInitialDisp[j]; 01074 } 01075 01076 if (nodeJInitialDisp != 0) { 01077 for (int j=0; j<3; j++) 01078 ug[j+3] -= nodeJInitialDisp[j]; 01079 } 01080 01081 // transform global end displacements to local coordinates 01082 static Vector ul(6); // total displacements 01083 01084 ul(0) = cosTheta*ug(0) + sinTheta*ug(1); 01085 ul(1) = -sinTheta*ug(0) + cosTheta*ug(1); 01086 ul(2) = ug(2); 01087 ul(3) = cosTheta*ug(3) + sinTheta*ug(4); 01088 ul(4) = -sinTheta*ug(3) + cosTheta*ug(4); 01089 ul(5) = ug(5); 01090 01091 if (nodeIOffset != 0) { 01092 double t02 = -cosTheta*nodeIOffset[1] + sinTheta*nodeIOffset[0]; 01093 double t12 = sinTheta*nodeIOffset[1] + cosTheta*nodeIOffset[0]; 01094 01095 ul(0) += t02*ug(2); 01096 ul(1) += t12*ug(2); 01097 } 01098 01099 if (nodeJOffset != 0) { 01100 double t35 = -cosTheta*nodeJOffset[1] + sinTheta*nodeJOffset[0]; 01101 double t45 = sinTheta*nodeJOffset[1] + cosTheta*nodeJOffset[0]; 01102 01103 ul(3) += t35*ug(5); 01104 ul(4) += t45*ug(5); 01105 } 01106 01107 // compute displacements at point xi, in local coordinates 01108 static Vector uxl(2), uxg(2); 01109 01110 uxl(0) = uxb(0) + ul(0); 01111 uxl(1) = uxb(1) + (1-xi)*ul(1) + xi*ul(4); 01112 01113 // rotate displacements to global coordinates 01114 // uxg = RljT*uxl 01115 uxg(0) = cosTheta*uxl(0) - sinTheta*uxl(1); 01116 uxg(1) = sinTheta*uxl(0) + cosTheta*uxl(1); 01117 01118 return uxg; 01119 } 01120 01121 01122 void 01123 LinearCrdTransf2d::Print(OPS_Stream &s, int flag) 01124 { 01125 s << "\nCrdTransf: " << this->getTag() << " Type: LinearCrdTransf2d"; 01126 if (nodeIOffset != 0) 01127 s << "\tnodeI Offset: " << nodeIOffset[0] << ' ' << nodeIOffset[1] << endln; 01128 if (nodeJOffset != 0) 01129 s << "\tnodeJ Offset: " << nodeJOffset[0] << ' ' << nodeJOffset[1] << endln; 01130 01131 } 01132 01133 01134 // AddingSensitivity:BEGIN /////////////////////////////// 01135 const Vector & 01136 LinearCrdTransf2d::getGlobalResistingForceShapeSensitivity(const Vector &pb, 01137 const Vector &p0, 01138 int gradNumber) 01139 { 01140 // transform resisting forces from the basic system to local coordinates 01141 static double pl[6]; 01142 01143 double q0 = pb(0); 01144 double q1 = pb(1); 01145 double q2 = pb(2); 01146 01147 double oneOverL = 1.0/L; 01148 01149 double V = oneOverL*(q1+q2); 01150 pl[0] = -q0; 01151 pl[1] = V; 01152 pl[2] = q1; 01153 pl[3] = q0; 01154 pl[4] = -V; 01155 pl[5] = q2; 01156 01157 // add end forces due to element p0 loads 01158 pl[0] += p0(0); 01159 pl[1] += p0(1); 01160 pl[4] += p0(2); 01161 01162 // transform resisting forces from local to global coordinates 01163 static Vector pg(6); 01164 pg.Zero(); 01165 01166 static ID nodeParameterID(2); 01167 nodeParameterID(0) = nodeIPtr->getCrdsSensitivity(); 01168 nodeParameterID(1) = nodeJPtr->getCrdsSensitivity(); 01169 01170 if (nodeParameterID(0) != 0 || nodeParameterID(1) != 0) { 01171 01172 if (nodeIOffset != 0 || nodeJOffset != 0) { 01173 opserr << "ERROR: Currently a node offset cannot be used in " << endln 01174 << " conjunction with random nodal coordinates." << endln; 01175 } 01176 01177 double dcosdh, dsindh, d1oLdh; 01178 01179 double dx = cosTheta*L; 01180 double dy = sinTheta*L; 01181 01182 if (nodeParameterID(0) == 1) { // here x1 is random 01183 dcosdh = (-L+dx*dx/L)/(L*L); 01184 dsindh = dx*dy/(L*L*L); 01185 d1oLdh = dx/(L*L*L); 01186 } 01187 if (nodeParameterID(0) == 2) { // here y1 is random 01188 dsindh = (-L+dy*dy/L)/(L*L); 01189 dcosdh = dx*dy/(L*L*L); 01190 d1oLdh = dy/(L*L*L); 01191 } 01192 01193 if (nodeParameterID(1) == 1) { // here x2 is random 01194 dcosdh = (L-dx*dx/L)/(L*L); 01195 dsindh = -dx*dy/(L*L*L); 01196 d1oLdh = -dx/(L*L*L); 01197 } 01198 if (nodeParameterID(1) == 2) { // here y2 is random 01199 dsindh = (L-dy*dy/L)/(L*L); 01200 dcosdh = -dx*dy/(L*L*L); 01201 d1oLdh = -dy/(L*L*L); 01202 } 01203 01204 pg(0) = dcosdh*pl[0] - dsindh*pl[1] - sinTheta*d1oLdh*(q1+q2); 01205 pg(1) = dsindh*pl[0] + dcosdh*pl[1] + cosTheta*d1oLdh*(q1+q2); 01206 01207 pg(3) = dcosdh*pl[3] - dsindh*pl[4] + sinTheta*d1oLdh*(q1+q2); 01208 pg(4) = dsindh*pl[3] + dcosdh*pl[4] - cosTheta*d1oLdh*(q1+q2); 01209 01210 pg(2) = 0.0; 01211 pg(5) = 0.0; 01212 } 01213 01214 return pg; 01215 } 01216 01217 const Vector & 01218 LinearCrdTransf2d::getBasicDisplSensitivity(int gradNumber) 01219 { 01220 static Vector U(6); 01221 static Vector dUdh(6); 01222 01223 const Vector &dispI = nodeIPtr->getTrialDisp(); 01224 const Vector &dispJ = nodeJPtr->getTrialDisp(); 01225 01226 for (int i = 0; i < 3; i++) { 01227 U(i) = dispI(i); 01228 U(i+3) = dispJ(i); 01229 dUdh(i) = nodeIPtr->getDispSensitivity((i+1),gradNumber); 01230 dUdh(i+3) = nodeJPtr->getDispSensitivity((i+1),gradNumber); 01231 } 01232 01233 static Vector dvdh(3); 01234 01235 static Matrix A(6,6); 01236 01237 A(0,0) = A(3,3) = cosTheta; 01238 A(0,1) = A(3,4) = sinTheta; 01239 A(1,0) = A(4,3) = -sinTheta; 01240 A(1,1) = A(4,4) = cosTheta; 01241 A(2,2) = A(5,5) = 1.0; 01242 01243 double dcosThetadh = 0.0; 01244 double dsinThetadh = 0.0; 01245 01246 double dx = cosTheta*L; 01247 double dy = sinTheta*L; 01248 01249 int nodeIid = nodeIPtr->getCrdsSensitivity(); 01250 int nodeJid = nodeJPtr->getCrdsSensitivity(); 01251 01252 //if (nodeIid == 0 && nodeJid == 0) 01253 // return dvdh; 01254 01255 if (nodeIid == 1) { // here x1 is random 01256 dcosThetadh = (-L+dx*dx/L)/(L*L); 01257 dsinThetadh = dx*dy/(L*L*L); 01258 } 01259 if (nodeIid == 2) { // here y1 is random 01260 dsinThetadh = (-L+dy*dy/L)/(L*L); 01261 dcosThetadh = dx*dy/(L*L*L); 01262 } 01263 01264 if (nodeJid == 1) { // here x2 is random 01265 dcosThetadh = (L-dx*dx/L)/(L*L); 01266 dsinThetadh = -dx*dy/(L*L*L); 01267 } 01268 if (nodeJid == 2) { // here y2 is random 01269 dsinThetadh = (L-dy*dy/L)/(L*L); 01270 dcosThetadh = -dx*dy/(L*L*L); 01271 } 01272 01273 static Matrix dAdh(6,6); 01274 01275 dAdh(0,0) = dAdh(3,3) = dcosThetadh; 01276 dAdh(1,1) = dAdh(4,4) = dcosThetadh; 01277 dAdh(0,1) = dAdh(3,4) = dsinThetadh; 01278 dAdh(1,0) = dAdh(4,3) = -dsinThetadh; 01279 01280 static Vector dudh(6); 01281 dudh = A*dUdh + dAdh*U; 01282 01283 static Vector u(6); 01284 u = A*U; 01285 01286 static Matrix Abl(3,6); 01287 01288 Abl(0,0) = -1.0; 01289 01290 Abl(1,1) = 1.0/L; 01291 Abl(2,1) = 1.0/L; 01292 01293 Abl(1,2) = 1; 01294 01295 Abl(0,3) = 1; 01296 01297 Abl(1,4) = -1.0/L; 01298 Abl(2,4) = -1.0/L; 01299 01300 Abl(2,5) = 1; 01301 01302 static Matrix dAbldh(3,6); 01303 01304 double dLdh = this->getdLdh(); 01305 double doneOverLdh = -dLdh/(L*L); 01306 01307 dAbldh(1,1) = doneOverLdh; 01308 dAbldh(2,1) = doneOverLdh; 01309 dAbldh(1,4) = -doneOverLdh; 01310 dAbldh(1,4) = -doneOverLdh; 01311 01312 dvdh = Abl*dudh + dAbldh*u; 01313 01314 return dvdh; 01315 } 01316 01317 bool 01318 LinearCrdTransf2d::isShapeSensitivity(void) 01319 { 01320 int nodeParameterI, nodeParameterJ; 01321 nodeParameterI = nodeIPtr->getCrdsSensitivity(); 01322 nodeParameterJ = nodeJPtr->getCrdsSensitivity(); 01323 01324 return (nodeParameterI != 0 || nodeParameterJ != 0); 01325 } 01326 01327 double 01328 LinearCrdTransf2d::getdLdh(void) 01329 { 01330 int nodeParameterI, nodeParameterJ; 01331 nodeParameterI = nodeIPtr->getCrdsSensitivity(); 01332 nodeParameterJ = nodeJPtr->getCrdsSensitivity(); 01333 01334 if (nodeParameterI != 0 || nodeParameterJ != 0) { 01335 01336 if (nodeIOffset != 0 || nodeJOffset != 0) { 01337 opserr << "ERROR: Currently a node offset cannot be used in " << endln 01338 << " conjunction with random nodal coordinates." << endln; 01339 } 01340 01341 if (nodeParameterI == 1) // here x1 is random 01342 return -cosTheta; 01343 if (nodeParameterI == 2) // here y1 is random 01344 return -sinTheta; 01345 01346 if (nodeParameterJ == 1) // here x2 is random 01347 return cosTheta; 01348 if (nodeParameterJ == 2) // here y2 is random 01349 return sinTheta; 01350 } 01351 else 01352 return 0.0; 01353 } 01354 01355 double 01356 LinearCrdTransf2d::getd1overLdh(void) 01357 { 01358 int nodeParameterI, nodeParameterJ; 01359 nodeParameterI = nodeIPtr->getCrdsSensitivity(); 01360 nodeParameterJ = nodeJPtr->getCrdsSensitivity(); 01361 01362 if (nodeParameterI != 0 || nodeParameterJ != 0) { 01363 01364 if (nodeIOffset != 0 || nodeJOffset != 0) { 01365 opserr << "ERROR: Currently a node offset cannot be used in " << endln 01366 << " conjunction with random nodal coordinates." << endln; 01367 } 01368 01369 if (nodeParameterI == 1) // here x1 is random 01370 return cosTheta/(L*L); 01371 if (nodeParameterI == 2) // here y1 is random 01372 return sinTheta/(L*L); 01373 01374 if (nodeParameterJ == 1) // here x2 is random 01375 return -cosTheta/(L*L); 01376 if (nodeParameterJ == 2) // here y2 is random 01377 return -sinTheta/(L*L); 01378 } 01379 else 01380 return 0.0; 01381 } 01382 // AddingSensitivity:END /////////////////////////////////////
Generated on Mon Oct 23 15:04:59 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents