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