Rev 1894 | Rev 2891 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 272 | mhscott | 1 | /* ****************************************************************** ** |
| 2 | ** OpenSees - Open System for Earthquake Engineering Simulation ** |
||
| 3 | ** Pacific Earthquake Engineering Research Center ** |
||
| 4 | ** ** |
||
| 5 | ** ** |
||
| 6 | ** (C) Copyright 1999, The Regents of the University of California ** |
||
| 7 | ** All Rights Reserved. ** |
||
| 8 | ** ** |
||
| 9 | ** Commercial use of this program without express permission of the ** |
||
| 10 | ** University of California, Berkeley, is strictly prohibited. See ** |
||
| 11 | ** file 'COPYRIGHT' in main directory for information on usage and ** |
||
| 12 | ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
||
| 13 | ** ** |
||
| 14 | ** Developed by: ** |
||
| 15 | ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
||
| 16 | ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
||
| 17 | ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
||
| 18 | ** ** |
||
| 19 | ** ****************************************************************** */ |
||
| 2221 | fmk | 20 | |
| 21 | // $Revision: 1.13 $ |
||
| 22 | // $Date: 2005-12-15 00:30:38 $ |
||
| 272 | mhscott | 23 | // $Source: /usr/local/cvs/OpenSees/SRC/coordTransformation/PDeltaCrdTransf3d.cpp,v $ |
| 2221 | fmk | 24 | |
| 25 | |||
| 272 | mhscott | 26 | // Written: Remo Magalhaes de Souza (rmsouza@ce.berkeley.edu) |
| 27 | // Created: 04/2000 |
||
| 28 | // Revision: A |
||
| 2221 | fmk | 29 | // |
| 30 | // Modified: 04/2005 Andreas Schellenberg (getBasicTrialVel, getBasicTrialAccel) |
||
| 272 | mhscott | 31 | // |
| 32 | // Purpose: This file contains the implementation for the |
||
| 33 | // PDeltaCrdTransf3d class. PDeltaCrdTransf3d is a linear |
||
| 34 | // transformation for a planar frame between the global |
||
| 35 | // and basic coordinate systems |
||
| 36 | |||
| 37 | |||
| 38 | #include <Vector.h> |
||
| 39 | #include <Matrix.h> |
||
| 40 | #include <Node.h> |
||
| 41 | #include <Channel.h> |
||
| 42 | |||
| 43 | #include <PDeltaCrdTransf3d.h> |
||
| 44 | |||
| 1843 | fmk | 45 | |
| 272 | mhscott | 46 | // constructor: |
| 47 | PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane): |
||
| 2221 | fmk | 48 | CrdTransf3d(tag, CRDTR_TAG_PDeltaCrdTransf3d), |
| 49 | nodeIPtr(0), nodeJPtr(0), |
||
| 50 | nodeIOffset(0), nodeJOffset(0), |
||
| 51 | L(0), ul17(0), ul28(0), |
||
| 52 | nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false) |
||
| 272 | mhscott | 53 | { |
| 2221 | fmk | 54 | for (int i = 0; i < 2; i++) |
| 55 | for (int j = 0; j < 3; j++) |
||
| 56 | R[i][j] = 0.0; |
||
| 57 | |||
| 58 | R[2][0] = vecInLocXZPlane(0); |
||
| 59 | R[2][1] = vecInLocXZPlane(1); |
||
| 60 | R[2][2] = vecInLocXZPlane(2); |
||
| 61 | |||
| 62 | // Does nothing |
||
| 63 | } |
||
| 272 | mhscott | 64 | |
| 610 | mhscott | 65 | |
| 272 | mhscott | 66 | // constructor: |
| 67 | PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane, |
||
| 2221 | fmk | 68 | const Vector &rigJntOffset1, |
| 69 | const Vector &rigJntOffset2): |
||
| 70 | CrdTransf3d(tag, CRDTR_TAG_PDeltaCrdTransf3d), |
||
| 71 | nodeIPtr(0), nodeJPtr(0), |
||
| 72 | nodeIOffset(0), nodeJOffset(0), |
||
| 73 | L(0), ul17(0), ul28(0), |
||
| 74 | nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false) |
||
| 272 | mhscott | 75 | { |
| 2221 | fmk | 76 | for (int i = 0; i < 2; i++) |
| 77 | for (int j = 0; j < 3; j++) |
||
| 78 | R[i][j] = 0.0; |
||
| 79 | |||
| 80 | R[2][0] = vecInLocXZPlane(0); |
||
| 81 | R[2][1] = vecInLocXZPlane(1); |
||
| 82 | R[2][2] = vecInLocXZPlane(2); |
||
| 83 | |||
| 84 | // check rigid joint offset for node I |
||
| 85 | if (&rigJntOffset1 == 0 || rigJntOffset1.Size() != 3 ) { |
||
| 86 | opserr << "PDeltaCrdTransf3d::PDeltaCrdTransf3d: Invalid rigid joint offset vector for node I\n"; |
||
| 87 | opserr << "Size must be 3\n"; |
||
| 88 | } |
||
| 89 | else if (rigJntOffset1.Norm() > 0.0) { |
||
| 90 | nodeIOffset = new double[3]; |
||
| 91 | nodeIOffset[0] = rigJntOffset1(0); |
||
| 92 | nodeIOffset[1] = rigJntOffset1(1); |
||
| 93 | nodeIOffset[2] = rigJntOffset1(2); |
||
| 94 | } |
||
| 95 | |||
| 96 | // check rigid joint offset for node J |
||
| 97 | if (&rigJntOffset2 == 0 || rigJntOffset2.Size() != 3 ) { |
||
| 98 | opserr << "PDeltaCrdTransf3d::PDeltaCrdTransf3d: Invalid rigid joint offset vector for node J\n"; |
||
| 99 | opserr << "Size must be 3\n"; |
||
| 100 | } |
||
| 101 | else if (rigJntOffset2.Norm() > 0.0) { |
||
| 102 | nodeJOffset = new double[3]; |
||
| 103 | nodeJOffset[0] = rigJntOffset2(0); |
||
| 104 | nodeJOffset[1] = rigJntOffset2(1); |
||
| 105 | nodeJOffset[2] = rigJntOffset2(2); |
||
| 106 | } |
||
| 272 | mhscott | 107 | } |
| 108 | |||
| 109 | |||
| 110 | // constructor: |
||
| 111 | // invoked by a FEM_ObjectBroker, recvSelf() needs to be invoked on this object. |
||
| 112 | PDeltaCrdTransf3d::PDeltaCrdTransf3d(): |
||
| 2221 | fmk | 113 | CrdTransf3d(0, CRDTR_TAG_PDeltaCrdTransf3d), |
| 114 | nodeIPtr(0), nodeJPtr(0), |
||
| 115 | nodeIOffset(0), nodeJOffset(0), |
||
| 116 | L(0), ul17(0), ul28(0), |
||
| 117 | nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false) |
||
| 272 | mhscott | 118 | { |
| 2221 | fmk | 119 | for (int i = 0; i < 3; i++) |
| 120 | for (int j = 0; j < 3; j++) |
||
| 121 | R[i][j] = 0.0; |
||
| 272 | mhscott | 122 | } |
| 123 | |||
| 124 | |||
| 125 | // destructor: |
||
| 126 | PDeltaCrdTransf3d::~PDeltaCrdTransf3d() |
||
| 127 | { |
||
| 2221 | fmk | 128 | if (nodeIOffset) |
| 129 | delete [] nodeIOffset; |
||
| 130 | if (nodeJOffset) |
||
| 131 | delete [] nodeJOffset; |
||
| 132 | if (nodeIInitialDisp != 0) |
||
| 133 | delete [] nodeIInitialDisp; |
||
| 134 | if (nodeJInitialDisp != 0) |
||
| 135 | delete [] nodeJInitialDisp; |
||
| 272 | mhscott | 136 | } |
| 137 | |||
| 138 | |||
| 139 | int |
||
| 140 | PDeltaCrdTransf3d::commitState(void) |
||
| 141 | { |
||
| 2221 | fmk | 142 | return 0; |
| 272 | mhscott | 143 | } |
| 144 | |||
| 145 | |||
| 146 | int |
||
| 147 | PDeltaCrdTransf3d::revertToLastCommit(void) |
||
| 148 | { |
||
| 2221 | fmk | 149 | return 0; |
| 272 | mhscott | 150 | } |
| 151 | |||
| 152 | |||
| 153 | int |
||
| 154 | PDeltaCrdTransf3d::revertToStart(void) |
||
| 155 | { |
||
| 2221 | fmk | 156 | return 0; |
| 272 | mhscott | 157 | } |
| 158 | |||
| 159 | |||
| 160 | int |
||
| 161 | PDeltaCrdTransf3d::initialize(Node *nodeIPointer, Node *nodeJPointer) |
||
| 162 | { |
||
| 2221 | fmk | 163 | int error; |
| 164 | |||
| 165 | nodeIPtr = nodeIPointer; |
||
| 166 | nodeJPtr = nodeJPointer; |
||
| 167 | |||
| 168 | if ((!nodeIPtr) || (!nodeJPtr)) |
||
| 169 | { |
||
| 170 | opserr << "\nPDeltaCrdTransf3d::initialize"; |
||
| 171 | opserr << "\ninvalid pointers to the element nodes\n"; |
||
| 172 | return -1; |
||
| 173 | } |
||
| 174 | |||
| 175 | // see if there is some initial displacements at nodes |
||
| 176 | if (initialDispChecked == false) { |
||
| 177 | const Vector &nodeIDisp = nodeIPtr->getDisp(); |
||
| 178 | const Vector &nodeJDisp = nodeJPtr->getDisp(); |
||
| 179 | for (int i=0; i<6; i++) |
||
| 180 | if (nodeIDisp(i) != 0.0) { |
||
| 181 | nodeIInitialDisp = new double [6]; |
||
| 182 | for (int j=0; j<6; j++) |
||
| 183 | nodeIInitialDisp[j] = nodeIDisp(j); |
||
| 184 | i = 6; |
||
| 185 | } |
||
| 186 | |||
| 187 | for (int j=0; j<6; j++) |
||
| 188 | if (nodeJDisp(j) != 0.0) { |
||
| 189 | nodeJInitialDisp = new double [6]; |
||
| 190 | for (int i=0; i<6; i++) |
||
| 191 | nodeJInitialDisp[i] = nodeJDisp(i); |
||
| 192 | j = 6; |
||
| 193 | } |
||
| 194 | |||
| 195 | initialDispChecked = true; |
||
| 196 | } |
||
| 197 | |||
| 198 | // get element length and orientation |
||
| 199 | if ((error = this->computeElemtLengthAndOrient())) |
||
| 200 | return error; |
||
| 201 | |||
| 202 | static Vector XAxis(3); |
||
| 203 | static Vector YAxis(3); |
||
| 204 | static Vector ZAxis(3); |
||
| 205 | |||
| 206 | // get 3by3 rotation matrix |
||
| 207 | if ((error = this->getLocalAxes(XAxis, YAxis, ZAxis))) |
||
| 208 | return error; |
||
| 209 | |||
| 210 | return 0; |
||
| 272 | mhscott | 211 | } |
| 212 | |||
| 213 | |||
| 214 | int |
||
| 215 | PDeltaCrdTransf3d::update(void) |
||
| 216 | { |
||
| 2221 | fmk | 217 | const Vector &disp1 = nodeIPtr->getTrialDisp(); |
| 218 | const Vector &disp2 = nodeJPtr->getTrialDisp(); |
||
| 1894 | fmk | 219 | |
| 2221 | fmk | 220 | static double ug[12]; |
| 221 | for (int i = 0; i < 6; i++) { |
||
| 222 | ug[i] = disp1(i); |
||
| 223 | ug[i+6] = disp2(i); |
||
| 224 | } |
||
| 1894 | fmk | 225 | |
| 2221 | fmk | 226 | if (nodeIInitialDisp != 0) { |
| 227 | for (int j=0; j<6; j++) |
||
| 228 | ug[j] -= nodeIInitialDisp[j]; |
||
| 229 | } |
||
| 230 | |||
| 231 | if (nodeJInitialDisp != 0) { |
||
| 232 | for (int j=0; j<6; j++) |
||
| 233 | ug[j+6] -= nodeJInitialDisp[j]; |
||
| 234 | } |
||
| 235 | |||
| 236 | double ul1, ul7, ul2, ul8; |
||
| 237 | |||
| 238 | ul1 = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2]; |
||
| 239 | ul2 = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2]; |
||
| 240 | |||
| 241 | ul7 = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8]; |
||
| 242 | ul8 = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8]; |
||
| 243 | |||
| 244 | static double Wu[3]; |
||
| 245 | |||
| 246 | if (nodeIOffset) { |
||
| 247 | Wu[0] = nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5]; |
||
| 248 | Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5]; |
||
| 249 | Wu[2] = nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4]; |
||
| 250 | |||
| 251 | ul1 += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 252 | ul2 += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 253 | } |
||
| 254 | |||
| 255 | if (nodeJOffset) { |
||
| 256 | Wu[0] = nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11]; |
||
| 257 | Wu[1] = -nodeJOffset[2]*ug[9] + nodeJOffset[0]*ug[11]; |
||
| 258 | Wu[2] = nodeJOffset[1]*ug[9] - nodeJOffset[0]*ug[10]; |
||
| 259 | |||
| 260 | ul7 += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 261 | ul8 += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 262 | } |
||
| 263 | |||
| 264 | ul17 = ul1-ul7; |
||
| 265 | ul28 = ul2-ul8; |
||
| 266 | |||
| 267 | return 0; |
||
| 272 | mhscott | 268 | } |
| 269 | |||
| 270 | |||
| 271 | int |
||
| 272 | PDeltaCrdTransf3d::computeElemtLengthAndOrient() |
||
| 273 | { |
||
| 2221 | fmk | 274 | // element projection |
| 275 | static Vector dx(3); |
||
| 276 | |||
| 277 | const Vector &ndICoords = nodeIPtr->getCrds(); |
||
| 278 | const Vector &ndJCoords = nodeJPtr->getCrds(); |
||
| 279 | |||
| 280 | dx(0) = ndJCoords(0) - ndICoords(0); |
||
| 281 | dx(1) = ndJCoords(1) - ndICoords(1); |
||
| 282 | dx(2) = ndJCoords(2) - ndICoords(2); |
||
| 283 | |||
| 284 | if (nodeIInitialDisp != 0) { |
||
| 285 | dx(0) -= nodeIInitialDisp[0]; |
||
| 286 | dx(1) -= nodeIInitialDisp[1]; |
||
| 287 | dx(2) -= nodeIInitialDisp[2]; |
||
| 288 | } |
||
| 289 | |||
| 290 | if (nodeJInitialDisp != 0) { |
||
| 291 | dx(0) += nodeJInitialDisp[0]; |
||
| 292 | dx(1) += nodeJInitialDisp[1]; |
||
| 293 | dx(2) += nodeJInitialDisp[2]; |
||
| 294 | } |
||
| 295 | |||
| 296 | if (nodeJOffset != 0) { |
||
| 297 | dx(0) += nodeJOffset[0]; |
||
| 298 | dx(1) += nodeJOffset[1]; |
||
| 299 | dx(2) += nodeJOffset[2]; |
||
| 300 | } |
||
| 301 | |||
| 302 | if (nodeIOffset != 0) { |
||
| 303 | dx(0) -= nodeIOffset[0]; |
||
| 304 | dx(1) -= nodeIOffset[1]; |
||
| 305 | dx(2) -= nodeIOffset[2]; |
||
| 306 | } |
||
| 307 | |||
| 308 | // calculate the element length |
||
| 309 | L = dx.Norm(); |
||
| 310 | |||
| 311 | if (L == 0.0) { |
||
| 312 | opserr << "\nPDeltaCrdTransf3d::computeElemtLengthAndOrien: 0 length\n"; |
||
| 313 | return -2; |
||
| 314 | } |
||
| 315 | |||
| 316 | // calculate the element local x axis components (direction cossines) |
||
| 317 | // wrt to the global coordinates |
||
| 318 | R[0][0] = dx(0)/L; |
||
| 319 | R[0][1] = dx(1)/L; |
||
| 320 | R[0][2] = dx(2)/L; |
||
| 321 | |||
| 322 | return 0; |
||
| 272 | mhscott | 323 | } |
| 324 | |||
| 325 | |||
| 326 | int |
||
| 1402 | fmk | 327 | PDeltaCrdTransf3d::getLocalAxes(Vector &XAxis, Vector &YAxis, Vector &ZAxis) |
| 272 | mhscott | 328 | { |
| 2221 | fmk | 329 | // Compute y = v cross x |
| 330 | // Note: v(i) is stored in R[2][i] |
||
| 331 | static Vector vAxis(3); |
||
| 332 | vAxis(0) = R[2][0]; vAxis(1) = R[2][1]; vAxis(2) = R[2][2]; |
||
| 333 | |||
| 334 | static Vector xAxis(3); |
||
| 335 | xAxis(0) = R[0][0]; xAxis(1) = R[0][1]; xAxis(2) = R[0][2]; |
||
| 336 | XAxis(0) = xAxis(0); XAxis(1) = xAxis(1); XAxis(2) = xAxis(2); |
||
| 337 | |||
| 338 | static Vector yAxis(3); |
||
| 339 | |||
| 340 | yAxis(0) = vAxis(1)*xAxis(2) - vAxis(2)*xAxis(1); |
||
| 341 | yAxis(1) = vAxis(2)*xAxis(0) - vAxis(0)*xAxis(2); |
||
| 342 | yAxis(2) = vAxis(0)*xAxis(1) - vAxis(1)*xAxis(0); |
||
| 343 | |||
| 344 | double ynorm = yAxis.Norm(); |
||
| 345 | |||
| 346 | if (ynorm == 0) { |
||
| 347 | opserr << "\nPDeltaCrdTransf3d::getLocalAxes"; |
||
| 348 | opserr << "\nvector v that defines plane xz is parallel to x axis\n"; |
||
| 349 | return -3; |
||
| 350 | } |
||
| 351 | |||
| 352 | yAxis /= ynorm; |
||
| 353 | |||
| 354 | YAxis(0) = yAxis(0); YAxis(1) = yAxis(1); YAxis(2) = yAxis(2); |
||
| 355 | |||
| 356 | // Compute z = x cross y |
||
| 357 | static Vector zAxis(3); |
||
| 358 | |||
| 359 | zAxis(0) = xAxis(1)*yAxis(2) - xAxis(2)*yAxis(1); |
||
| 360 | zAxis(1) = xAxis(2)*yAxis(0) - xAxis(0)*yAxis(2); |
||
| 361 | zAxis(2) = xAxis(0)*yAxis(1) - xAxis(1)*yAxis(0); |
||
| 362 | ZAxis(0) = zAxis(0); ZAxis(1) = zAxis(1); ZAxis(2) = zAxis(2); |
||
| 363 | |||
| 364 | // Fill in transformation matrix |
||
| 365 | R[1][0] = yAxis(0); |
||
| 366 | R[1][1] = yAxis(1); |
||
| 367 | R[1][2] = yAxis(2); |
||
| 368 | |||
| 369 | R[2][0] = zAxis(0); |
||
| 370 | R[2][1] = zAxis(1); |
||
| 371 | R[2][2] = zAxis(2); |
||
| 372 | |||
| 373 | return 0; |
||
| 272 | mhscott | 374 | } |
| 375 | |||
| 376 | |||
| 377 | double |
||
| 378 | PDeltaCrdTransf3d::getInitialLength(void) |
||
| 379 | { |
||
| 2221 | fmk | 380 | return L; |
| 272 | mhscott | 381 | } |
| 382 | |||
| 383 | |||
| 384 | double |
||
| 385 | PDeltaCrdTransf3d::getDeformedLength(void) |
||
| 386 | { |
||
| 2221 | fmk | 387 | return L; |
| 272 | mhscott | 388 | } |
| 389 | |||
| 390 | |||
| 391 | const Vector & |
||
| 392 | PDeltaCrdTransf3d::getBasicTrialDisp (void) |
||
| 393 | { |
||
| 2221 | fmk | 394 | // determine global displacements |
| 395 | const Vector &disp1 = nodeIPtr->getTrialDisp(); |
||
| 396 | const Vector &disp2 = nodeJPtr->getTrialDisp(); |
||
| 1894 | fmk | 397 | |
| 2221 | fmk | 398 | static double ug[12]; |
| 399 | for (int i = 0; i < 6; i++) { |
||
| 400 | ug[i] = disp1(i); |
||
| 401 | ug[i+6] = disp2(i); |
||
| 402 | } |
||
| 1894 | fmk | 403 | |
| 2221 | fmk | 404 | if (nodeIInitialDisp != 0) { |
| 405 | for (int j=0; j<6; j++) |
||
| 406 | ug[j] -= nodeIInitialDisp[j]; |
||
| 407 | } |
||
| 408 | |||
| 409 | if (nodeJInitialDisp != 0) { |
||
| 410 | for (int j=0; j<6; j++) |
||
| 411 | ug[j+6] -= nodeJInitialDisp[j]; |
||
| 412 | } |
||
| 413 | |||
| 414 | double oneOverL = 1.0/L; |
||
| 415 | |||
| 416 | static Vector ub(6); |
||
| 417 | |||
| 418 | static double ul[12]; |
||
| 419 | |||
| 420 | ul[0] = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2]; |
||
| 421 | ul[1] = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2]; |
||
| 422 | ul[2] = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2]; |
||
| 423 | |||
| 424 | ul[3] = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5]; |
||
| 425 | ul[4] = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5]; |
||
| 426 | ul[5] = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5]; |
||
| 427 | |||
| 428 | ul[6] = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8]; |
||
| 429 | ul[7] = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8]; |
||
| 430 | ul[8] = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8]; |
||
| 431 | |||
| 432 | ul[9] = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11]; |
||
| 433 | ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11]; |
||
| 434 | ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11]; |
||
| 435 | |||
| 436 | static double Wu[3]; |
||
| 437 | if (nodeIOffset) { |
||
| 438 | Wu[0] = nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5]; |
||
| 439 | Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5]; |
||
| 440 | Wu[2] = nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4]; |
||
| 441 | |||
| 442 | ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 443 | ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 444 | ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 445 | } |
||
| 446 | |||
| 447 | if (nodeJOffset) { |
||
| 448 | Wu[0] = nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11]; |
||
| 449 | Wu[1] = -nodeJOffset[2]*ug[9] + nodeJOffset[0]*ug[11]; |
||
| 450 | Wu[2] = nodeJOffset[1]*ug[9] - nodeJOffset[0]*ug[10]; |
||
| 451 | |||
| 452 | ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 453 | ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 454 | ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 455 | } |
||
| 456 | |||
| 457 | ub(0) = ul[6] - ul[0]; |
||
| 458 | double tmp; |
||
| 459 | tmp = oneOverL*(ul[1]-ul[7]); |
||
| 460 | ub(1) = ul[5] + tmp; |
||
| 461 | ub(2) = ul[11] + tmp; |
||
| 462 | tmp = oneOverL*(ul[8]-ul[2]); |
||
| 463 | ub(3) = ul[4] + tmp; |
||
| 464 | ub(4) = ul[10] + tmp; |
||
| 465 | ub(5) = ul[9] - ul[3]; |
||
| 466 | |||
| 467 | return ub; |
||
| 272 | mhscott | 468 | } |
| 469 | |||
| 470 | |||
| 471 | const Vector & |
||
| 472 | PDeltaCrdTransf3d::getBasicIncrDisp (void) |
||
| 473 | { |
||
| 2221 | fmk | 474 | // determine global displacements |
| 475 | const Vector &disp1 = nodeIPtr->getIncrDisp(); |
||
| 476 | const Vector &disp2 = nodeJPtr->getIncrDisp(); |
||
| 477 | |||
| 478 | static double ug[12]; |
||
| 479 | for (int i = 0; i < 6; i++) { |
||
| 480 | ug[i] = disp1(i); |
||
| 481 | ug[i+6] = disp2(i); |
||
| 482 | } |
||
| 483 | |||
| 484 | double oneOverL = 1.0/L; |
||
| 485 | |||
| 486 | static Vector ub(6); |
||
| 487 | |||
| 488 | static double ul[12]; |
||
| 489 | |||
| 490 | ul[0] = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2]; |
||
| 491 | ul[1] = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2]; |
||
| 492 | ul[2] = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2]; |
||
| 493 | |||
| 494 | ul[3] = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5]; |
||
| 495 | ul[4] = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5]; |
||
| 496 | ul[5] = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5]; |
||
| 497 | |||
| 498 | ul[6] = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8]; |
||
| 499 | ul[7] = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8]; |
||
| 500 | ul[8] = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8]; |
||
| 501 | |||
| 502 | ul[9] = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11]; |
||
| 503 | ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11]; |
||
| 504 | ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11]; |
||
| 505 | |||
| 506 | static double Wu[3]; |
||
| 507 | if (nodeIOffset) { |
||
| 508 | Wu[0] = nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5]; |
||
| 509 | Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5]; |
||
| 510 | Wu[2] = nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4]; |
||
| 511 | |||
| 512 | ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 513 | ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 514 | ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 515 | } |
||
| 516 | |||
| 517 | if (nodeJOffset) { |
||
| 518 | Wu[0] = nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11]; |
||
| 519 | Wu[1] = -nodeJOffset[2]*ug[9] + nodeJOffset[0]*ug[11]; |
||
| 520 | Wu[2] = nodeJOffset[1]*ug[9] - nodeJOffset[0]*ug[10]; |
||
| 521 | |||
| 522 | ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 523 | ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 524 | ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 525 | } |
||
| 526 | |||
| 527 | ub(0) = ul[6] - ul[0]; |
||
| 528 | double tmp; |
||
| 529 | tmp = oneOverL*(ul[1]-ul[7]); |
||
| 530 | ub(1) = ul[5] + tmp; |
||
| 531 | ub(2) = ul[11] + tmp; |
||
| 532 | tmp = oneOverL*(ul[8]-ul[2]); |
||
| 533 | ub(3) = ul[4] + tmp; |
||
| 534 | ub(4) = ul[10] + tmp; |
||
| 535 | ub(5) = ul[9] - ul[3]; |
||
| 536 | |||
| 537 | return ub; |
||
| 538 | } |
||
| 272 | mhscott | 539 | |
| 2221 | fmk | 540 | |
| 541 | const Vector & |
||
| 542 | PDeltaCrdTransf3d::getBasicIncrDeltaDisp(void) |
||
| 543 | { |
||
| 544 | // determine global displacements |
||
| 545 | const Vector &disp1 = nodeIPtr->getIncrDeltaDisp(); |
||
| 546 | const Vector &disp2 = nodeJPtr->getIncrDeltaDisp(); |
||
| 547 | |||
| 548 | static double ug[12]; |
||
| 549 | for (int i = 0; i < 6; i++) { |
||
| 550 | ug[i] = disp1(i); |
||
| 551 | ug[i+6] = disp2(i); |
||
| 552 | } |
||
| 553 | |||
| 554 | double oneOverL = 1.0/L; |
||
| 555 | |||
| 556 | static Vector ub(6); |
||
| 557 | |||
| 558 | static double ul[12]; |
||
| 559 | |||
| 560 | ul[0] = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2]; |
||
| 561 | ul[1] = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2]; |
||
| 562 | ul[2] = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2]; |
||
| 563 | |||
| 564 | ul[3] = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5]; |
||
| 565 | ul[4] = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5]; |
||
| 566 | ul[5] = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5]; |
||
| 567 | |||
| 568 | ul[6] = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8]; |
||
| 569 | ul[7] = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8]; |
||
| 570 | ul[8] = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8]; |
||
| 571 | |||
| 572 | ul[9] = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11]; |
||
| 573 | ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11]; |
||
| 574 | ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11]; |
||
| 575 | |||
| 576 | static double Wu[3]; |
||
| 577 | if (nodeIOffset) { |
||
| 578 | Wu[0] = nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5]; |
||
| 579 | Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5]; |
||
| 580 | Wu[2] = nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4]; |
||
| 581 | |||
| 582 | ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 583 | ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 584 | ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 585 | } |
||
| 586 | |||
| 587 | if (nodeJOffset) { |
||
| 588 | Wu[0] = nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11]; |
||
| 589 | Wu[1] = -nodeJOffset[2]*ug[9] + nodeJOffset[0]*ug[11]; |
||
| 590 | Wu[2] = nodeJOffset[1]*ug[9] - nodeJOffset[0]*ug[10]; |
||
| 591 | |||
| 592 | ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 593 | ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 594 | ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 595 | } |
||
| 596 | |||
| 597 | ub(0) = ul[6] - ul[0]; |
||
| 598 | double tmp; |
||
| 599 | tmp = oneOverL*(ul[1]-ul[7]); |
||
| 600 | ub(1) = ul[5] + tmp; |
||
| 601 | ub(2) = ul[11] + tmp; |
||
| 602 | tmp = oneOverL*(ul[8]-ul[2]); |
||
| 603 | ub(3) = ul[4] + tmp; |
||
| 604 | ub(4) = ul[10] + tmp; |
||
| 605 | ub(5) = ul[9] - ul[3]; |
||
| 606 | |||
| 607 | return ub; |
||
| 608 | } |
||
| 609 | |||
| 610 | |||
| 611 | const Vector & |
||
| 612 | PDeltaCrdTransf3d::getBasicTrialVel(void) |
||
| 613 | { |
||
| 614 | // determine global velocities |
||
| 615 | const Vector &vel1 = nodeIPtr->getTrialVel(); |
||
| 616 | const Vector &vel2 = nodeJPtr->getTrialVel(); |
||
| 617 | |||
| 618 | static double vg[12]; |
||
| 272 | mhscott | 619 | for (int i = 0; i < 6; i++) { |
| 2221 | fmk | 620 | vg[i] = vel1(i); |
| 621 | vg[i+6] = vel2(i); |
||
| 272 | mhscott | 622 | } |
| 2221 | fmk | 623 | |
| 272 | mhscott | 624 | double oneOverL = 1.0/L; |
| 2221 | fmk | 625 | |
| 626 | static Vector vb(6); |
||
| 627 | |||
| 628 | static double vl[12]; |
||
| 629 | |||
| 630 | vl[0] = R[0][0]*vg[0] + R[0][1]*vg[1] + R[0][2]*vg[2]; |
||
| 631 | vl[1] = R[1][0]*vg[0] + R[1][1]*vg[1] + R[1][2]*vg[2]; |
||
| 632 | vl[2] = R[2][0]*vg[0] + R[2][1]*vg[1] + R[2][2]*vg[2]; |
||
| 633 | |||
| 634 | vl[3] = R[0][0]*vg[3] + R[0][1]*vg[4] + R[0][2]*vg[5]; |
||
| 635 | vl[4] = R[1][0]*vg[3] + R[1][1]*vg[4] + R[1][2]*vg[5]; |
||
| 636 | vl[5] = R[2][0]*vg[3] + R[2][1]*vg[4] + R[2][2]*vg[5]; |
||
| 637 | |||
| 638 | vl[6] = R[0][0]*vg[6] + R[0][1]*vg[7] + R[0][2]*vg[8]; |
||
| 639 | vl[7] = R[1][0]*vg[6] + R[1][1]*vg[7] + R[1][2]*vg[8]; |
||
| 640 | vl[8] = R[2][0]*vg[6] + R[2][1]*vg[7] + R[2][2]*vg[8]; |
||
| 641 | |||
| 642 | vl[9] = R[0][0]*vg[9] + R[0][1]*vg[10] + R[0][2]*vg[11]; |
||
| 643 | vl[10] = R[1][0]*vg[9] + R[1][1]*vg[10] + R[1][2]*vg[11]; |
||
| 644 | vl[11] = R[2][0]*vg[9] + R[2][1]*vg[10] + R[2][2]*vg[11]; |
||
| 645 | |||
| 272 | mhscott | 646 | static double Wu[3]; |
| 647 | if (nodeIOffset) { |
||
| 2221 | fmk | 648 | Wu[0] = nodeIOffset[2]*vg[4] - nodeIOffset[1]*vg[5]; |
| 649 | Wu[1] = -nodeIOffset[2]*vg[3] + nodeIOffset[0]*vg[5]; |
||
| 650 | Wu[2] = nodeIOffset[1]*vg[3] - nodeIOffset[0]*vg[4]; |
||
| 651 | |||
| 652 | vl[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 653 | vl[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 654 | vl[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 272 | mhscott | 655 | } |
| 2221 | fmk | 656 | |
| 272 | mhscott | 657 | if (nodeJOffset) { |
| 2221 | fmk | 658 | Wu[0] = nodeJOffset[2]*vg[10] - nodeJOffset[1]*vg[11]; |
| 659 | Wu[1] = -nodeJOffset[2]*vg[9] + nodeJOffset[0]*vg[11]; |
||
| 660 | Wu[2] = nodeJOffset[1]*vg[9] - nodeJOffset[0]*vg[10]; |
||
| 661 | |||
| 662 | vl[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 663 | vl[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 664 | vl[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 272 | mhscott | 665 | } |
| 2221 | fmk | 666 | |
| 667 | vb(0) = vl[6] - vl[0]; |
||
| 272 | mhscott | 668 | double tmp; |
| 2221 | fmk | 669 | tmp = oneOverL*(vl[1]-vl[7]); |
| 670 | vb(1) = vl[5] + tmp; |
||
| 671 | vb(2) = vl[11] + tmp; |
||
| 672 | tmp = oneOverL*(vl[8]-vl[2]); |
||
| 673 | vb(3) = vl[4] + tmp; |
||
| 674 | vb(4) = vl[10] + tmp; |
||
| 675 | vb(5) = vl[9] - vl[3]; |
||
| 676 | |||
| 677 | return vb; |
||
| 272 | mhscott | 678 | } |
| 679 | |||
| 680 | |||
| 681 | const Vector & |
||
| 2221 | fmk | 682 | PDeltaCrdTransf3d::getBasicTrialAccel(void) |
| 272 | mhscott | 683 | { |
| 2221 | fmk | 684 | // determine global accelerations |
| 685 | const Vector &accel1 = nodeIPtr->getTrialAccel(); |
||
| 686 | const Vector &accel2 = nodeJPtr->getTrialAccel(); |
||
| 687 | |||
| 688 | static double ag[12]; |
||
| 272 | mhscott | 689 | for (int i = 0; i < 6; i++) { |
| 2221 | fmk | 690 | ag[i] = accel1(i); |
| 691 | ag[i+6] = accel2(i); |
||
| 272 | mhscott | 692 | } |
| 2221 | fmk | 693 | |
| 272 | mhscott | 694 | double oneOverL = 1.0/L; |
| 2221 | fmk | 695 | |
| 696 | static Vector ab(6); |
||
| 697 | |||
| 698 | static double al[12]; |
||
| 699 | |||
| 700 | al[0] = R[0][0]*ag[0] + R[0][1]*ag[1] + R[0][2]*ag[2]; |
||
| 701 | al[1] = R[1][0]*ag[0] + R[1][1]*ag[1] + R[1][2]*ag[2]; |
||
| 702 | al[2] = R[2][0]*ag[0] + R[2][1]*ag[1] + R[2][2]*ag[2]; |
||
| 703 | |||
| 704 | al[3] = R[0][0]*ag[3] + R[0][1]*ag[4] + R[0][2]*ag[5]; |
||
| 705 | al[4] = R[1][0]*ag[3] + R[1][1]*ag[4] + R[1][2]*ag[5]; |
||
| 706 | al[5] = R[2][0]*ag[3] + R[2][1]*ag[4] + R[2][2]*ag[5]; |
||
| 707 | |||
| 708 | al[6] = R[0][0]*ag[6] + R[0][1]*ag[7] + R[0][2]*ag[8]; |
||
| 709 | al[7] = R[1][0]*ag[6] + R[1][1]*ag[7] + R[1][2]*ag[8]; |
||
| 710 | al[8] = R[2][0]*ag[6] + R[2][1]*ag[7] + R[2][2]*ag[8]; |
||
| 711 | |||
| 712 | al[9] = R[0][0]*ag[9] + R[0][1]*ag[10] + R[0][2]*ag[11]; |
||
| 713 | al[10] = R[1][0]*ag[9] + R[1][1]*ag[10] + R[1][2]*ag[11]; |
||
| 714 | al[11] = R[2][0]*ag[9] + R[2][1]*ag[10] + R[2][2]*ag[11]; |
||
| 715 | |||
| 272 | mhscott | 716 | static double Wu[3]; |
| 717 | if (nodeIOffset) { |
||
| 2221 | fmk | 718 | Wu[0] = nodeIOffset[2]*ag[4] - nodeIOffset[1]*ag[5]; |
| 719 | Wu[1] = -nodeIOffset[2]*ag[3] + nodeIOffset[0]*ag[5]; |
||
| 720 | Wu[2] = nodeIOffset[1]*ag[3] - nodeIOffset[0]*ag[4]; |
||
| 721 | |||
| 722 | al[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 723 | al[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 724 | al[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 272 | mhscott | 725 | } |
| 2221 | fmk | 726 | |
| 272 | mhscott | 727 | if (nodeJOffset) { |
| 2221 | fmk | 728 | Wu[0] = nodeJOffset[2]*ag[10] - nodeJOffset[1]*ag[11]; |
| 729 | Wu[1] = -nodeJOffset[2]*ag[9] + nodeJOffset[0]*ag[11]; |
||
| 730 | Wu[2] = nodeJOffset[1]*ag[9] - nodeJOffset[0]*ag[10]; |
||
| 731 | |||
| 732 | al[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 733 | al[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 734 | al[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 272 | mhscott | 735 | } |
| 2221 | fmk | 736 | |
| 737 | ab(0) = al[6] - al[0]; |
||
| 272 | mhscott | 738 | double tmp; |
| 2221 | fmk | 739 | tmp = oneOverL*(al[1]-al[7]); |
| 740 | ab(1) = al[5] + tmp; |
||
| 741 | ab(2) = al[11] + tmp; |
||
| 742 | tmp = oneOverL*(al[8]-al[2]); |
||
| 743 | ab(3) = al[4] + tmp; |
||
| 744 | ab(4) = al[10] + tmp; |
||
| 745 | ab(5) = al[9] - al[3]; |
||
| 746 | |||
| 747 | return ab; |
||
| 272 | mhscott | 748 | } |
| 749 | |||
| 750 | |||
| 751 | const Vector & |
||
| 1002 | mhscott | 752 | PDeltaCrdTransf3d::getGlobalResistingForce(const Vector &pb, const Vector &p0) |
| 272 | mhscott | 753 | { |
| 2221 | fmk | 754 | // transform resisting forces from the basic system to local coordinates |
| 755 | static double pl[12]; |
||
| 756 | |||
| 757 | double q0 = pb(0); |
||
| 758 | double q1 = pb(1); |
||
| 759 | double q2 = pb(2); |
||
| 760 | double q3 = pb(3); |
||
| 761 | double q4 = pb(4); |
||
| 762 | double q5 = pb(5); |
||
| 763 | |||
| 764 | double oneOverL = 1.0/L; |
||
| 765 | |||
| 766 | pl[0] = -q0; |
||
| 767 | pl[1] = oneOverL*(q1+q2); |
||
| 768 | pl[2] = -oneOverL*(q3+q4); |
||
| 769 | pl[3] = -q5; |
||
| 770 | pl[4] = q3; |
||
| 771 | pl[5] = q1; |
||
| 772 | pl[6] = q0; |
||
| 773 | pl[7] = -pl[1]; |
||
| 774 | pl[8] = -pl[2]; |
||
| 775 | pl[9] = q5; |
||
| 776 | pl[10] = q4; |
||
| 777 | pl[11] = q2; |
||
| 778 | |||
| 779 | pl[0] += p0(0); |
||
| 780 | pl[1] += p0(1); |
||
| 781 | pl[7] += p0(2); |
||
| 782 | pl[2] += p0(3); |
||
| 783 | pl[8] += p0(4); |
||
| 784 | |||
| 785 | // Include leaning column effects (P-Delta) |
||
| 786 | double NoverL; |
||
| 787 | NoverL = ul17*q0*oneOverL; |
||
| 788 | pl[1] += NoverL; |
||
| 789 | pl[7] -= NoverL; |
||
| 790 | NoverL = ul28*q0*oneOverL; |
||
| 791 | pl[2] += NoverL; |
||
| 792 | pl[8] -= NoverL; |
||
| 793 | |||
| 794 | // transform resisting forces from local to global coordinates |
||
| 795 | static Vector pg(12); |
||
| 796 | |||
| 797 | pg(0) = R[0][0]*pl[0] + R[1][0]*pl[1] + R[2][0]*pl[2]; |
||
| 798 | pg(1) = R[0][1]*pl[0] + R[1][1]*pl[1] + R[2][1]*pl[2]; |
||
| 799 | pg(2) = R[0][2]*pl[0] + R[1][2]*pl[1] + R[2][2]*pl[2]; |
||
| 800 | |||
| 801 | pg(3) = R[0][0]*pl[3] + R[1][0]*pl[4] + R[2][0]*pl[5]; |
||
| 802 | pg(4) = R[0][1]*pl[3] + R[1][1]*pl[4] + R[2][1]*pl[5]; |
||
| 803 | pg(5) = R[0][2]*pl[3] + R[1][2]*pl[4] + R[2][2]*pl[5]; |
||
| 804 | |||
| 805 | pg(6) = R[0][0]*pl[6] + R[1][0]*pl[7] + R[2][0]*pl[8]; |
||
| 806 | pg(7) = R[0][1]*pl[6] + R[1][1]*pl[7] + R[2][1]*pl[8]; |
||
| 807 | pg(8) = R[0][2]*pl[6] + R[1][2]*pl[7] + R[2][2]*pl[8]; |
||
| 808 | |||
| 809 | pg(9) = R[0][0]*pl[9] + R[1][0]*pl[10] + R[2][0]*pl[11]; |
||
| 810 | pg(10) = R[0][1]*pl[9] + R[1][1]*pl[10] + R[2][1]*pl[11]; |
||
| 811 | pg(11) = R[0][2]*pl[9] + R[1][2]*pl[10] + R[2][2]*pl[11]; |
||
| 812 | |||
| 813 | if (nodeIOffset) { |
||
| 814 | pg(3) += -nodeIOffset[2]*pg(1) + nodeIOffset[1]*pg(2); |
||
| 815 | pg(4) += nodeIOffset[2]*pg(0) - nodeIOffset[0]*pg(2); |
||
| 816 | pg(5) += -nodeIOffset[1]*pg(0) + nodeIOffset[0]*pg(1); |
||
| 817 | } |
||
| 818 | |||
| 819 | if (nodeJOffset) { |
||
| 820 | pg(9) += -nodeJOffset[2]*pg(7) + nodeJOffset[1]*pg(8); |
||
| 821 | pg(10) += nodeJOffset[2]*pg(6) - nodeJOffset[0]*pg(8); |
||
| 822 | pg(11) += -nodeJOffset[1]*pg(6) + nodeJOffset[0]*pg(7); |
||
| 823 | } |
||
| 824 | |||
| 825 | return pg; |
||
| 272 | mhscott | 826 | } |
| 827 | |||
| 1843 | fmk | 828 | |
| 272 | mhscott | 829 | const Matrix & |
| 830 | PDeltaCrdTransf3d::getGlobalStiffMatrix (const Matrix &KB, const Vector &pb) |
||
| 831 | { |
||
| 2221 | fmk | 832 | static Matrix kg(12,12); // Global stiffness for return |
| 833 | static double kb[6][6]; // Basic stiffness |
||
| 834 | static double kl[12][12]; // Local stiffness |
||
| 835 | static double tmp[12][12]; // Temporary storage |
||
| 836 | |||
| 837 | double oneOverL = 1.0/L; |
||
| 838 | |||
| 839 | int i,j; |
||
| 840 | for (i = 0; i < 6; i++) |
||
| 841 | for (j = 0; j < 6; j++) |
||
| 842 | kb[i][j] = KB(i,j); |
||
| 843 | |||
| 844 | // Transform basic stiffness to local system |
||
| 845 | // First compute kb*T_{bl} |
||
| 846 | for (i = 0; i < 6; i++) { |
||
| 847 | tmp[i][0] = -kb[i][0]; |
||
| 848 | tmp[i][1] = oneOverL*(kb[i][1]+kb[i][2]); |
||
| 849 | tmp[i][2] = -oneOverL*(kb[i][3]+kb[i][4]); |
||
| 850 | tmp[i][3] = -kb[i][5]; |
||
| 851 | tmp[i][4] = kb[i][3]; |
||
| 852 | tmp[i][5] = kb[i][1]; |
||
| 853 | tmp[i][6] = kb[i][0]; |
||
| 854 | tmp[i][7] = -tmp[i][1]; |
||
| 855 | tmp[i][8] = -tmp[i][2]; |
||
| 856 | tmp[i][9] = kb[i][5]; |
||
| 857 | tmp[i][10] = kb[i][4]; |
||
| 858 | tmp[i][11] = kb[i][2]; |
||
| 859 | } |
||
| 860 | |||
| 861 | // Now compute T'_{bl}*(kb*T_{bl}) |
||
| 862 | for (i = 0; i < 12; i++) { |
||
| 863 | kl[0][i] = -tmp[0][i]; |
||
| 864 | kl[1][i] = oneOverL*(tmp[1][i]+tmp[2][i]); |
||
| 865 | kl[2][i] = -oneOverL*(tmp[3][i]+tmp[4][i]); |
||
| 866 | kl[3][i] = -tmp[5][i]; |
||
| 867 | kl[4][i] = tmp[3][i]; |
||
| 868 | kl[5][i] = tmp[1][i]; |
||
| 869 | kl[6][i] = tmp[0][i]; |
||
| 870 | kl[7][i] = -kl[1][i]; |
||
| 871 | kl[8][i] = -kl[2][i]; |
||
| 872 | kl[9][i] = tmp[5][i]; |
||
| 873 | kl[10][i] = tmp[4][i]; |
||
| 874 | kl[11][i] = tmp[2][i]; |
||
| 875 | } |
||
| 876 | |||
| 877 | // Include geometric stiffness effects in local system |
||
| 878 | double NoverL = pb(0)*oneOverL; |
||
| 879 | kl[1][1] += NoverL; |
||
| 880 | kl[2][2] += NoverL; |
||
| 881 | kl[7][7] += NoverL; |
||
| 882 | kl[8][8] += NoverL; |
||
| 883 | kl[1][7] -= NoverL; |
||
| 884 | kl[7][1] -= NoverL; |
||
| 885 | kl[2][8] -= NoverL; |
||
| 886 | kl[8][2] -= NoverL; |
||
| 887 | |||
| 888 | static double RWI[3][3]; |
||
| 889 | |||
| 890 | if (nodeIOffset) { |
||
| 891 | // Compute RWI |
||
| 892 | RWI[0][0] = -R[0][1]*nodeIOffset[2] + R[0][2]*nodeIOffset[1]; |
||
| 893 | RWI[1][0] = -R[1][1]*nodeIOffset[2] + R[1][2]*nodeIOffset[1]; |
||
| 894 | RWI[2][0] = -R[2][1]*nodeIOffset[2] + R[2][2]*nodeIOffset[1]; |
||
| 895 | |||
| 896 | RWI[0][1] = R[0][0]*nodeIOffset[2] - R[0][2]*nodeIOffset[0]; |
||
| 897 | RWI[1][1] = R[1][0]*nodeIOffset[2] - R[1][2]*nodeIOffset[0]; |
||
| 898 | RWI[2][1] = R[2][0]*nodeIOffset[2] - R[2][2]*nodeIOffset[0]; |
||
| 899 | |||
| 900 | RWI[0][2] = -R[0][0]*nodeIOffset[1] + R[0][1]*nodeIOffset[0]; |
||
| 901 | RWI[1][2] = -R[1][0]*nodeIOffset[1] + R[1][1]*nodeIOffset[0]; |
||
| 902 | RWI[2][2] = -R[2][0]*nodeIOffset[1] + R[2][1]*nodeIOffset[0]; |
||
| 903 | } |
||
| 904 | |||
| 905 | static double RWJ[3][3]; |
||
| 906 | |||
| 907 | if (nodeJOffset) { |
||
| 908 | // Compute RWJ |
||
| 909 | RWJ[0][0] = -R[0][1]*nodeJOffset[2] + R[0][2]*nodeJOffset[1]; |
||
| 910 | RWJ[1][0] = -R[1][1]*nodeJOffset[2] + R[1][2]*nodeJOffset[1]; |
||
| 911 | RWJ[2][0] = -R[2][1]*nodeJOffset[2] + R[2][2]*nodeJOffset[1]; |
||
| 912 | |||
| 913 | RWJ[0][1] = R[0][0]*nodeJOffset[2] - R[0][2]*nodeJOffset[0]; |
||
| 914 | RWJ[1][1] = R[1][0]*nodeJOffset[2] - R[1][2]*nodeJOffset[0]; |
||
| 915 | RWJ[2][1] = R[2][0]*nodeJOffset[2] - R[2][2]*nodeJOffset[0]; |
||
| 916 | |||
| 917 | RWJ[0][2] = -R[0][0]*nodeJOffset[1] + R[0][1]*nodeJOffset[0]; |
||
| 918 | RWJ[1][2] = -R[1][0]*nodeJOffset[1] + R[1][1]*nodeJOffset[0]; |
||
| 919 | RWJ[2][2] = -R[2][0]*nodeJOffset[1] + R[2][1]*nodeJOffset[0]; |
||
| 920 | } |
||
| 921 | |||
| 922 | // Transform local stiffness to global system |
||
| 923 | // First compute kl*T_{lg} |
||
| 924 | int m; |
||
| 925 | for (m = 0; m < 12; m++) { |
||
| 926 | tmp[m][0] = kl[m][0]*R[0][0] + kl[m][1]*R[1][0] + kl[m][2]*R[2][0]; |
||
| 927 | tmp[m][1] = kl[m][0]*R[0][1] + kl[m][1]*R[1][1] + kl[m][2]*R[2][1]; |
||
| 928 | tmp[m][2] = kl[m][0]*R[0][2] + kl[m][1]*R[1][2] + kl[m][2]*R[2][2]; |
||
| 929 | |||
| 930 | tmp[m][3] = kl[m][3]*R[0][0] + kl[m][4]*R[1][0] + kl[m][5]*R[2][0]; |
||
| 931 | tmp[m][4] = kl[m][3]*R[0][1] + kl[m][4]*R[1][1] + kl[m][5]*R[2][1]; |
||
| 932 | tmp[m][5] = kl[m][3]*R[0][2] + kl[m][4]*R[1][2] + kl[m][5]*R[2][2]; |
||
| 933 | |||
| 934 | if (nodeIOffset) { |
||
| 935 | tmp[m][3] += kl[m][0]*RWI[0][0] + kl[m][1]*RWI[1][0] + kl[m][2]*RWI[2][0]; |
||
| 936 | tmp[m][4] += kl[m][0]*RWI[0][1] + kl[m][1]*RWI[1][1] + kl[m][2]*RWI[2][1]; |
||
| 937 | tmp[m][5] += kl[m][0]*RWI[0][2] + kl[m][1]*RWI[1][2] + kl[m][2]*RWI[2][2]; |
||
| 938 | } |
||
| 939 | |||
| 940 | tmp[m][6] = kl[m][6]*R[0][0] + kl[m][7]*R[1][0] + kl[m][8]*R[2][0]; |
||
| 941 | tmp[m][7] = kl[m][6]*R[0][1] + kl[m][7]*R[1][1] + kl[m][8]*R[2][1]; |
||
| 942 | tmp[m][8] = kl[m][6]*R[0][2] + kl[m][7]*R[1][2] + kl[m][8]*R[2][2]; |
||
| 943 | |||
| 944 | tmp[m][9] = kl[m][9]*R[0][0] + kl[m][10]*R[1][0] + kl[m][11]*R[2][0]; |
||
| 945 | tmp[m][10] = kl[m][9]*R[0][1] + kl[m][10]*R[1][1] + kl[m][11]*R[2][1]; |
||
| 946 | tmp[m][11] = kl[m][9]*R[0][2] + kl[m][10]*R[1][2] + kl[m][11]*R[2][2]; |
||
| 947 | |||
| 948 | if (nodeJOffset) { |
||
| 949 | tmp[m][9] += kl[m][6]*RWJ[0][0] + kl[m][7]*RWJ[1][0] + kl[m][8]*RWJ[2][0]; |
||
| 950 | tmp[m][10] += kl[m][6]*RWJ[0][1] + kl[m][7]*RWJ[1][1] + kl[m][8]*RWJ[2][1]; |
||
| 951 | tmp[m][11] += kl[m][6]*RWJ[0][2] + kl[m][7]*RWJ[1][2] + kl[m][8]*RWJ[2][2]; |
||
| 952 | } |
||
| 953 | |||
| 954 | } |
||
| 955 | |||
| 956 | // Now compute T'_{lg}*(kl*T_{lg}) |
||
| 957 | for (m = 0; m < 12; m++) { |
||
| 958 | kg(0,m) = R[0][0]*tmp[0][m] + R[1][0]*tmp[1][m] + R[2][0]*tmp[2][m]; |
||
| 959 | kg(1,m) = R[0][1]*tmp[0][m] + R[1][1]*tmp[1][m] + R[2][1]*tmp[2][m]; |
||
| 960 | kg(2,m) = R[0][2]*tmp[0][m] + R[1][2]*tmp[1][m] + R[2][2]*tmp[2][m]; |
||
| 961 | |||
| 962 | kg(3,m) = R[0][0]*tmp[3][m] + R[1][0]*tmp[4][m] + R[2][0]*tmp[5][m]; |
||
| 963 | kg(4,m) = R[0][1]*tmp[3][m] + R[1][1]*tmp[4][m] + R[2][1]*tmp[5][m]; |
||
| 964 | kg(5,m) = R[0][2]*tmp[3][m] + R[1][2]*tmp[4][m] + R[2][2]*tmp[5][m]; |
||
| 965 | |||
| 966 | if (nodeIOffset) { |
||
| 967 | kg(3,m) += RWI[0][0]*tmp[0][m] + RWI[1][0]*tmp[1][m] + RWI[2][0]*tmp[2][m]; |
||
| 968 | kg(4,m) += RWI[0][1]*tmp[0][m] + RWI[1][1]*tmp[1][m] + RWI[2][1]*tmp[2][m]; |
||
| 969 | kg(5,m) += RWI[0][2]*tmp[0][m] + RWI[1][2]*tmp[1][m] + RWI[2][2]*tmp[2][m]; |
||
| 970 | } |
||
| 971 | |||
| 972 | kg(6,m) = R[0][0]*tmp[6][m] + R[1][0]*tmp[7][m] + R[2][0]*tmp[8][m]; |
||
| 973 | kg(7,m) = R[0][1]*tmp[6][m] + R[1][1]*tmp[7][m] + R[2][1]*tmp[8][m]; |
||
| 974 | kg(8,m) = R[0][2]*tmp[6][m] + R[1][2]*tmp[7][m] + R[2][2]*tmp[8][m]; |
||
| 975 | |||
| 976 | kg(9,m) = R[0][0]*tmp[9][m] + R[1][0]*tmp[10][m] + R[2][0]*tmp[11][m]; |
||
| 977 | kg(10,m) = R[0][1]*tmp[9][m] + R[1][1]*tmp[10][m] + R[2][1]*tmp[11][m]; |
||
| 978 | kg(11,m) = R[0][2]*tmp[9][m] + R[1][2]*tmp[10][m] + R[2][2]*tmp[11][m]; |
||
| 979 | |||
| 980 | if (nodeJOffset) { |
||
| 981 | kg(9,m) += RWJ[0][0]*tmp[6][m] + RWJ[1][0]*tmp[7][m] + RWJ[2][0]*tmp[8][m]; |
||
| 982 | kg(10,m) += RWJ[0][1]*tmp[6][m] + RWJ[1][1]*tmp[7][m] + RWJ[2][1]*tmp[8][m]; |
||
| 983 | kg(11,m) += RWJ[0][2]*tmp[6][m] + RWJ[1][2]*tmp[7][m] + RWJ[2][2]*tmp[8][m]; |
||
| 984 | } |
||
| 985 | } |
||
| 986 | |||
| 987 | return kg; |
||
| 272 | mhscott | 988 | } |
| 1108 | fmk | 989 | |
| 990 | |||
| 991 | const Matrix & |
||
| 992 | PDeltaCrdTransf3d::getInitialGlobalStiffMatrix (const Matrix &KB) |
||
| 993 | { |
||
| 2221 | fmk | 994 | static Matrix kg(12,12); // Global stiffness for return |
| 995 | static double kb[6][6]; // Basic stiffness |
||
| 996 | static double kl[12][12]; // Local stiffness |
||
| 997 | static double tmp[12][12]; // Temporary storage |
||
| 998 | |||
| 999 | double oneOverL = 1.0/L; |
||
| 1000 | |||
| 1001 | int i,j; |
||
| 1002 | for (i = 0; i < 6; i++) |
||
| 1003 | for (j = 0; j < 6; j++) |
||
| 1004 | kb[i][j] = KB(i,j); |
||
| 1005 | |||
| 1006 | // Transform basic stiffness to local system |
||
| 1007 | // First compute kb*T_{bl} |
||
| 1008 | for (i = 0; i < 6; i++) { |
||
| 1009 | tmp[i][0] = -kb[i][0]; |
||
| 1010 | tmp[i][1] = oneOverL*(kb[i][1]+kb[i][2]); |
||
| 1011 | tmp[i][2] = -oneOverL*(kb[i][3]+kb[i][4]); |
||
| 1012 | tmp[i][3] = -kb[i][5]; |
||
| 1013 | tmp[i][4] = kb[i][3]; |
||
| 1014 | tmp[i][5] = kb[i][1]; |
||
| 1015 | tmp[i][6] = kb[i][0]; |
||
| 1016 | tmp[i][7] = -tmp[i][1]; |
||
| 1017 | tmp[i][8] = -tmp[i][2]; |
||
| 1018 | tmp[i][9] = kb[i][5]; |
||
| 1019 | tmp[i][10] = kb[i][4]; |
||
| 1020 | tmp[i][11] = kb[i][2]; |
||
| 1021 | } |
||
| 1022 | |||
| 1023 | // Now compute T'_{bl}*(kb*T_{bl}) |
||
| 1024 | for (i = 0; i < 12; i++) { |
||
| 1025 | kl[0][i] = -tmp[0][i]; |
||
| 1026 | kl[1][i] = oneOverL*(tmp[1][i]+tmp[2][i]); |
||
| 1027 | kl[2][i] = -oneOverL*(tmp[3][i]+tmp[4][i]); |
||
| 1028 | kl[3][i] = -tmp[5][i]; |
||
| 1029 | kl[4][i] = tmp[3][i]; |
||
| 1030 | kl[5][i] = tmp[1][i]; |
||
| 1031 | kl[6][i] = tmp[0][i]; |
||
| 1032 | kl[7][i] = -kl[1][i]; |
||
| 1033 | kl[8][i] = -kl[2][i]; |
||
| 1034 | kl[9][i] = tmp[5][i]; |
||
| 1035 | kl[10][i] = tmp[4][i]; |
||
| 1036 | kl[11][i] = tmp[2][i]; |
||
| 1037 | } |
||
| 1038 | |||
| 1039 | // Include geometric stiffness effects in local system |
||
| 1040 | // double NoverL = pb(0)*oneOverL; |
||
| 1041 | //kl[1][1] += NoverL; |
||
| 1042 | //kl[2][2] += NoverL; |
||
| 1043 | //kl[7][7] += NoverL; |
||
| 1044 | //kl[8][8] += NoverL; |
||
| 1045 | //kl[1][7] -= NoverL; |
||
| 1046 | //kl[7][1] -= NoverL; |
||
| 1047 | //kl[2][8] -= NoverL; |
||
| 1048 | //kl[8][2] -= NoverL; |
||
| 1049 | |||
| 1050 | |||
| 1051 | static double RWI[3][3]; |
||
| 1052 | |||
| 1053 | if (nodeIOffset) { |
||
| 1054 | // Compute RWI |
||
| 1055 | RWI[0][0] = -R[0][1]*nodeIOffset[2] + R[0][2]*nodeIOffset[1]; |
||
| 1056 | RWI[1][0] = -R[1][1]*nodeIOffset[2] + R[1][2]*nodeIOffset[1]; |
||
| 1057 | RWI[2][0] = -R[2][1]*nodeIOffset[2] + R[2][2]*nodeIOffset[1]; |
||
| 1058 | |||
| 1059 | RWI[0][1] = R[0][0]*nodeIOffset[2] - R[0][2]*nodeIOffset[0]; |
||
| 1060 | RWI[1][1] = R[1][0]*nodeIOffset[2] - R[1][2]*nodeIOffset[0]; |
||
| 1061 | RWI[2][1] = R[2][0]*nodeIOffset[2] - R[2][2]*nodeIOffset[0]; |
||
| 1062 | |||
| 1063 | RWI[0][2] = -R[0][0]*nodeIOffset[1] + R[0][1]*nodeIOffset[0]; |
||
| 1064 | RWI[1][2] = -R[1][0]*nodeIOffset[1] + R[1][1]*nodeIOffset[0]; |
||
| 1065 | RWI[2][2] = -R[2][0]*nodeIOffset[1] + R[2][1]*nodeIOffset[0]; |
||
| 1066 | } |
||
| 1067 | |||
| 1068 | static double RWJ[3][3]; |
||
| 1069 | |||
| 1070 | if (nodeJOffset) { |
||
| 1071 | // Compute RWJ |
||
| 1072 | RWJ[0][0] = -R[0][1]*nodeJOffset[2] + R[0][2]*nodeJOffset[1]; |
||
| 1073 | RWJ[1][0] = -R[1][1]*nodeJOffset[2] + R[1][2]*nodeJOffset[1]; |
||
| 1074 | RWJ[2][0] = -R[2][1]*nodeJOffset[2] + R[2][2]*nodeJOffset[1]; |
||
| 1075 | |||
| 1076 | RWJ[0][1] = R[0][0]*nodeJOffset[2] - R[0][2]*nodeJOffset[0]; |
||
| 1077 | RWJ[1][1] = R[1][0]*nodeJOffset[2] - R[1][2]*nodeJOffset[0]; |
||
| 1078 | RWJ[2][1] = R[2][0]*nodeJOffset[2] - R[2][2]*nodeJOffset[0]; |
||
| 1079 | |||
| 1080 | RWJ[0][2] = -R[0][0]*nodeJOffset[1] + R[0][1]*nodeJOffset[0]; |
||
| 1081 | RWJ[1][2] = -R[1][0]*nodeJOffset[1] + R[1][1]*nodeJOffset[0]; |
||
| 1082 | RWJ[2][2] = -R[2][0]*nodeJOffset[1] + R[2][1]*nodeJOffset[0]; |
||
| 1083 | } |
||
| 1084 | |||
| 1085 | // Transform local stiffness to global system |
||
| 1086 | // First compute kl*T_{lg} |
||
| 1087 | |||
| 1088 | int m; |
||
| 1089 | for (m = 0; m < 12; m++) { |
||
| 1090 | tmp[m][0] = kl[m][0]*R[0][0] + kl[m][1]*R[1][0] + kl[m][2]*R[2][0]; |
||
| 1091 | tmp[m][1] = kl[m][0]*R[0][1] + kl[m][1]*R[1][1] + kl[m][2]*R[2][1]; |
||
| 1092 | tmp[m][2] = kl[m][0]*R[0][2] + kl[m][1]*R[1][2] + kl[m][2]*R[2][2]; |
||
| 1093 | |||
| 1094 | tmp[m][3] = kl[m][3]*R[0][0] + kl[m][4]*R[1][0] + kl[m][5]*R[2][0]; |
||
| 1095 | tmp[m][4] = kl[m][3]*R[0][1] + kl[m][4]*R[1][1] + kl[m][5]*R[2][1]; |
||
| 1096 | tmp[m][5] = kl[m][3]*R[0][2] + kl[m][4]*R[1][2] + kl[m][5]*R[2][2]; |
||
| 1097 | |||
| 1098 | if (nodeIOffset) { |
||
| 1099 | tmp[m][3] += kl[m][0]*RWI[0][0] + kl[m][1]*RWI[1][0] + kl[m][2]*RWI[2][0]; |
||
| 1100 | tmp[m][4] += kl[m][0]*RWI[0][1] + kl[m][1]*RWI[1][1] + kl[m][2]*RWI[2][1]; |
||
| 1101 | tmp[m][5] += kl[m][0]*RWI[0][2] + kl[m][1]*RWI[1][2] + kl[m][2]*RWI[2][2]; |
||
| 1102 | } |
||
| 1103 | |||
| 1104 | tmp[m][6] = kl[m][6]*R[0][0] + kl[m][7]*R[1][0] + kl[m][8]*R[2][0]; |
||
| 1105 | tmp[m][7] = kl[m][6]*R[0][1] + kl[m][7]*R[1][1] + kl[m][8]*R[2][1]; |
||
| 1106 | tmp[m][8] = kl[m][6]*R[0][2] + kl[m][7]*R[1][2] + kl[m][8]*R[2][2]; |
||
| 1107 | |||
| 1108 | tmp[m][9] = kl[m][9]*R[0][0] + kl[m][10]*R[1][0] + kl[m][11]*R[2][0]; |
||
| 1109 | tmp[m][10] = kl[m][9]*R[0][1] + kl[m][10]*R[1][1] + kl[m][11]*R[2][1]; |
||
| 1110 | tmp[m][11] = kl[m][9]*R[0][2] + kl[m][10]*R[1][2] + kl[m][11]*R[2][2]; |
||
| 1111 | |||
| 1112 | if (nodeJOffset) { |
||
| 1113 | tmp[m][9] += kl[m][6]*RWJ[0][0] + kl[m][7]*RWJ[1][0] + kl[m][8]*RWJ[2][0]; |
||
| 1114 | tmp[m][10] += kl[m][6]*RWJ[0][1] + kl[m][7]*RWJ[1][1] + kl[m][8]*RWJ[2][1]; |
||
| 1115 | tmp[m][11] += kl[m][6]*RWJ[0][2] + kl[m][7]*RWJ[1][2] + kl[m][8]*RWJ[2][2]; |
||
| 1116 | } |
||
| 1117 | |||
| 1118 | } |
||
| 1119 | |||
| 1120 | // Now compute T'_{lg}*(kl*T_{lg}) |
||
| 1121 | for (m = 0; m < 12; m++) { |
||
| 1122 | kg(0,m) = R[0][0]*tmp[0][m] + R[1][0]*tmp[1][m] + R[2][0]*tmp[2][m]; |
||
| 1123 | kg(1,m) = R[0][1]*tmp[0][m] + R[1][1]*tmp[1][m] + R[2][1]*tmp[2][m]; |
||
| 1124 | kg(2,m) = R[0][2]*tmp[0][m] + R[1][2]*tmp[1][m] + R[2][2]*tmp[2][m]; |
||
| 1125 | |||
| 1126 | kg(3,m) = R[0][0]*tmp[3][m] + R[1][0]*tmp[4][m] + R[2][0]*tmp[5][m]; |
||
| 1127 | kg(4,m) = R[0][1]*tmp[3][m] + R[1][1]*tmp[4][m] + R[2][1]*tmp[5][m]; |
||
| 1128 | kg(5,m) = R[0][2]*tmp[3][m] + R[1][2]*tmp[4][m] + R[2][2]*tmp[5][m]; |
||
| 1129 | |||
| 1130 | if (nodeIOffset) { |
||
| 1131 | kg(3,m) += RWI[0][0]*tmp[0][m] + RWI[1][0]*tmp[1][m] + RWI[2][0]*tmp[2][m]; |
||
| 1132 | kg(4,m) += RWI[0][1]*tmp[0][m] + RWI[1][1]*tmp[1][m] + RWI[2][1]*tmp[2][m]; |
||
| 1133 | kg(5,m) += RWI[0][2]*tmp[0][m] + RWI[1][2]*tmp[1][m] + RWI[2][2]*tmp[2][m]; |
||
| 1134 | } |
||
| 1135 | |||
| 1136 | kg(6,m) = R[0][0]*tmp[6][m] + R[1][0]*tmp[7][m] + R[2][0]*tmp[8][m]; |
||
| 1137 | kg(7,m) = R[0][1]*tmp[6][m] + R[1][1]*tmp[7][m] + R[2][1]*tmp[8][m]; |
||
| 1138 | kg(8,m) = R[0][2]*tmp[6][m] + R[1][2]*tmp[7][m] + R[2][2]*tmp[8][m]; |
||
| 1139 | |||
| 1140 | kg(9,m) = R[0][0]*tmp[9][m] + R[1][0]*tmp[10][m] + R[2][0]*tmp[11][m]; |
||
| 1141 | kg(10,m) = R[0][1]*tmp[9][m] + R[1][1]*tmp[10][m] + R[2][1]*tmp[11][m]; |
||
| 1142 | kg(11,m) = R[0][2]*tmp[9][m] + R[1][2]*tmp[10][m] + R[2][2]*tmp[11][m]; |
||
| 1143 | |||
| 1144 | if (nodeJOffset) { |
||
| 1145 | kg(9,m) += RWJ[0][0]*tmp[6][m] + RWJ[1][0]*tmp[7][m] + RWJ[2][0]*tmp[8][m]; |
||
| 1146 | kg(10,m) += RWJ[0][1]*tmp[6][m] + RWJ[1][1]*tmp[7][m] + RWJ[2][1]*tmp[8][m]; |
||
| 1147 | kg(11,m) += RWJ[0][2]*tmp[6][m] + RWJ[1][2]*tmp[7][m] + RWJ[2][2]*tmp[8][m]; |
||
| 1148 | } |
||
| 1149 | } |
||
| 1150 | |||
| 1151 | return kg; |
||
| 1152 | } |
||
| 1108 | fmk | 1153 | |
| 1154 | |||
| 272 | mhscott | 1155 | CrdTransf3d * |
| 1156 | PDeltaCrdTransf3d::getCopy(void) |
||
| 1157 | { |
||
| 2221 | fmk | 1158 | // create a new instance of PDeltaCrdTransf3d |
| 1159 | |||
| 1160 | PDeltaCrdTransf3d *theCopy; |
||
| 1161 | |||
| 1162 | static Vector xz(3); |
||
| 1163 | xz(0) = R[2][0]; |
||
| 1164 | xz(1) = R[2][1]; |
||
| 1165 | xz(2) = R[2][2]; |
||
| 1166 | |||
| 1167 | Vector offsetI(3); |
||
| 1168 | Vector offsetJ(3); |
||
| 1169 | |||
| 1170 | if (nodeIOffset) { |
||
| 1171 | offsetI(0) = nodeIOffset[0]; |
||
| 1172 | offsetI(1) = nodeIOffset[1]; |
||
| 1173 | offsetI(2) = nodeIOffset[2]; |
||
| 1174 | } |
||
| 1175 | |||
| 1176 | if (nodeJOffset) { |
||
| 1177 | offsetJ(0) = nodeJOffset[0]; |
||
| 1178 | offsetJ(1) = nodeJOffset[1]; |
||
| 1179 | offsetJ(2) = nodeJOffset[2]; |
||
| 1180 | } |
||
| 1181 | |||
| 1182 | theCopy = new PDeltaCrdTransf3d(this->getTag(), xz, offsetI, offsetJ); |
||
| 1183 | |||
| 1184 | theCopy->nodeIPtr = nodeIPtr; |
||
| 1185 | theCopy->nodeJPtr = nodeJPtr; |
||
| 1186 | theCopy->L = L; |
||
| 1187 | theCopy->ul17 = ul17; |
||
| 1188 | theCopy->ul28 = ul28; |
||
| 1189 | for (int i = 0; i < 3; i++) |
||
| 1190 | for (int j = 0; j < 3; j++) |
||
| 1191 | theCopy->R[i][j] = R[i][j]; |
||
| 1192 | |||
| 1193 | |||
| 1194 | return theCopy; |
||
| 272 | mhscott | 1195 | } |
| 1196 | |||
| 1197 | |||
| 1198 | int |
||
| 1199 | PDeltaCrdTransf3d::sendSelf(int cTag, Channel &theChannel) |
||
| 1200 | { |
||
| 2221 | fmk | 1201 | int res = 0; |
| 1202 | |||
| 1203 | static Vector data(23); |
||
| 1204 | data(0) = this->getTag(); |
||
| 1205 | data(1) = L; |
||
| 1206 | |||
| 1207 | if (nodeIOffset != 0) { |
||
| 1208 | data(2) = nodeIOffset[0]; |
||
| 1209 | data(3) = nodeIOffset[1]; |
||
| 1210 | data(4) = nodeIOffset[2]; |
||
| 1211 | } else { |
||
| 1212 | data(2) = 0.0; |
||
| 1213 | data(3) = 0.0; |
||
| 1214 | data(4) = 0.0; |
||
| 1215 | } |
||
| 1216 | |||
| 1217 | if (nodeJOffset != 0) { |
||
| 1218 | data(5) = nodeJOffset[0]; |
||
| 1219 | data(6) = nodeJOffset[1]; |
||
| 1220 | data(7) = nodeJOffset[2]; |
||
| 1221 | } else { |
||
| 1222 | data(5) = 0.0; |
||
| 1223 | data(6) = 0.0; |
||
| 1224 | data(7) = 0.0; |
||
| 1225 | } |
||
| 1226 | |||
| 1227 | if (nodeIInitialDisp != 0) { |
||
| 1228 | data(8) = nodeIInitialDisp[0]; |
||
| 1229 | data(9) = nodeIInitialDisp[1]; |
||
| 1230 | data(10) = nodeIInitialDisp[2]; |
||
| 1231 | data(11) = nodeIInitialDisp[3]; |
||
| 1232 | data(12) = nodeIInitialDisp[4]; |
||
| 1233 | data(13) = nodeIInitialDisp[5]; |
||
| 1234 | } else { |
||
| 1235 | data(8) = 0.0; |
||
| 1236 | data(9) = 0.0; |
||
| 1237 | data(10) = 0.0; |
||
| 1238 | data(11) = 0.0; |
||
| 1239 | data(12) = 0.0; |
||
| 1240 | data(13) = 0.0; |
||
| 1241 | } |
||
| 1242 | |||
| 1243 | if (nodeJInitialDisp != 0) { |
||
| 1244 | data(14) = nodeJInitialDisp[0]; |
||
| 1245 | data(15) = nodeJInitialDisp[1]; |
||
| 1246 | data(16) = nodeJInitialDisp[2]; |
||
| 1247 | data(17) = nodeJInitialDisp[3]; |
||
| 1248 | data(18) = nodeJInitialDisp[4]; |
||
| 1249 | data(19) = nodeJInitialDisp[5]; |
||
| 1250 | } else { |
||
| 1251 | data(14) = 0.0; |
||
| 1252 | data(15) = 0.0; |
||
| 1253 | data(16) = 0.0; |
||
| 1254 | data(17) = 0.0; |
||
| 1255 | data(18) = 0.0; |
||
| 1256 | data(19) = 0.0; |
||
| 1257 | } |
||
| 1258 | |||
| 1259 | data(20) = R[2][0]; |
||
| 1260 | data(21) = R[2][1]; |
||
| 1261 | data(22) = R[2][2]; |
||
| 1262 | |||
| 1263 | res += theChannel.sendVector(this->getDbTag(), cTag, data); |
||
| 1264 | if (res < 0) { |
||
| 1265 | opserr << "PDeltaCrdTransf3d::sendSelf - failed to send Vector\n"; |
||
| 1266 | return res; |
||
| 1267 | } |
||
| 1268 | |||
| 272 | mhscott | 1269 | return res; |
| 1270 | } |
||
| 1271 | |||
| 1272 | |||
| 1273 | int |
||
| 1274 | PDeltaCrdTransf3d::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker) |
||
| 1275 | { |
||
| 2221 | fmk | 1276 | int res = 0; |
| 1277 | |||
| 1278 | static Vector data(13); |
||
| 1279 | |||
| 1280 | res += theChannel.recvVector(this->getDbTag(), cTag, data); |
||
| 1281 | if (res < 0) { |
||
| 1282 | opserr << "PDeltaCrdTransf3d::recvSelf - failed to receive Vector\n"; |
||
| 1283 | return res; |
||
| 1284 | } |
||
| 1285 | |||
| 1286 | this->setTag((int)data(0)); |
||
| 1287 | L = data(1); |
||
| 1288 | data(0) = this->getTag(); |
||
| 1289 | data(1) = L; |
||
| 1290 | |||
| 1291 | int flag; |
||
| 1292 | int i,j; |
||
| 1293 | |||
| 1294 | flag = 0; |
||
| 1295 | for (i=2; i<=4; i++) |
||
| 1296 | if (data(i) != 0.0) |
||
| 1297 | flag = 1; |
||
| 1298 | if (flag == 1) { |
||
| 1299 | if (nodeIOffset == 0) |
||
| 1300 | nodeIOffset = new double[3]; |
||
| 1301 | for (i=2, j=0; i<=4; i++, j++) |
||
| 1302 | nodeIOffset[j] = data(i); |
||
| 1303 | } |
||
| 1304 | |||
| 1305 | flag = 0; |
||
| 1306 | for (i=5; i<=7; i++) |
||
| 1307 | if (data(i) != 0.0) |
||
| 1308 | flag = 1; |
||
| 1309 | if (flag == 1) { |
||
| 1310 | if (nodeJOffset == 0) |
||
| 1311 | nodeJOffset = new double[3]; |
||
| 1312 | for (i=5, j=0; i<=7; i++, j++) |
||
| 1313 | nodeJOffset[j] = data(i); |
||
| 1314 | } |
||
| 1315 | |||
| 1316 | flag = 0; |
||
| 1317 | for (i=8; i<=13; i++) |
||
| 1318 | if (data(i) != 0.0) |
||
| 1319 | flag = 1; |
||
| 1320 | if (flag == 1) { |
||
| 1321 | if (nodeIInitialDisp == 0) |
||
| 1322 | nodeIInitialDisp = new double[6]; |
||
| 1323 | for (i=8, j=0; i<=13; i++, j++) |
||
| 1324 | nodeIInitialDisp[j] = data(i); |
||
| 1325 | } |
||
| 1326 | |||
| 1327 | flag = 0; |
||
| 1328 | for (i=14; i<=19; i++) |
||
| 1329 | if (data(i) != 0.0) |
||
| 1330 | flag = 1; |
||
| 1331 | if (flag == 1) { |
||
| 1332 | if (nodeJInitialDisp == 0) |
||
| 1333 | nodeJInitialDisp = new double [6]; |
||
| 1334 | for (i=14, j=0; i<=19; i++, j++) |
||
| 1335 | nodeJInitialDisp[j] = data(i); |
||
| 1336 | } |
||
| 1337 | |||
| 1338 | R[2][0] = data(20); |
||
| 1339 | R[2][1] = data(21); |
||
| 1340 | R[2][2] = data(22); |
||
| 1341 | |||
| 1342 | initialDispChecked = true; |
||
| 1343 | |||
| 1344 | return res; |
||
| 1345 | } |
||
| 1894 | fmk | 1346 | |
| 1347 | |||
| 272 | mhscott | 1348 | const Vector & |
| 1349 | PDeltaCrdTransf3d::getPointGlobalCoordFromLocal(const Vector &xl) |
||
| 1350 | { |
||
| 2221 | fmk | 1351 | static Vector xg(3); |
| 1352 | |||
| 1353 | //xg = nodeIPtr->getCrds() + nodeIOffset; |
||
| 1354 | xg = nodeIPtr->getCrds(); |
||
| 1355 | |||
| 1356 | if (nodeIOffset) { |
||
| 1357 | xg(0) += nodeIOffset[0]; |
||
| 1358 | xg(1) += nodeIOffset[1]; |
||
| 1359 | xg(2) += nodeIOffset[2]; |
||
| 1360 | } |
||
| 1361 | |||
| 1362 | // xg = xg + Rlj'*xl |
||
| 1363 | //xg.addMatrixTransposeVector(1.0, Rlj, xl, 1.0); |
||
| 1364 | xg(0) += R[0][0]*xl(0) + R[1][0]*xl(1) + R[2][0]*xl(2); |
||
| 1365 | xg(1) += R[0][1]*xl(0) + R[1][1]*xl(1) + R[2][1]*xl(2); |
||
| 1366 | xg(2) += R[0][2]*xl(0) + R[1][2]*xl(1) + R[2][2]*xl(2); |
||
| 1367 | |||
| 1368 | return xg; |
||
| 1369 | } |
||
| 272 | mhscott | 1370 | |
| 1371 | |||
| 1372 | const Vector & |
||
| 1373 | PDeltaCrdTransf3d::getPointGlobalDisplFromBasic (double xi, const Vector &uxb) |
||
| 1374 | { |
||
| 2221 | fmk | 1375 | // determine global displacements |
| 1376 | const Vector &disp1 = nodeIPtr->getTrialDisp(); |
||
| 1377 | const Vector &disp2 = nodeJPtr->getTrialDisp(); |
||
| 1378 | |||
| 1379 | static double ug[12]; |
||
| 1380 | for (int i = 0; i < 6; i++) |
||
| 1381 | { |
||
| 1382 | ug[i] = disp1(i); |
||
| 1383 | ug[i+6] = disp2(i); |
||
| 1384 | } |
||
| 1385 | |||
| 1386 | if (nodeIInitialDisp != 0) { |
||
| 1387 | for (int j=0; j<6; j++) |
||
| 1388 | ug[j] -= nodeIInitialDisp[j]; |
||
| 1389 | } |
||
| 1390 | |||
| 1391 | if (nodeJInitialDisp != 0) { |
||
| 1392 | for (int j=0; j<6; j++) |
||
| 1393 | ug[j+6] -= nodeJInitialDisp[j]; |
||
| 1394 | } |
||
| 1395 | |||
| 1396 | // transform global end displacements to local coordinates |
||
| 1397 | //ul.addMatrixVector(0.0, Tlg, ug, 1.0); // ul = Tlg * ug; |
||
| 1398 | static double ul[12]; |
||
| 1399 | |||
| 1400 | ul[0] = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2]; |
||
| 1401 | ul[1] = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2]; |
||
| 1402 | ul[2] = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2]; |
||
| 1403 | |||
| 1404 | ul[7] = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8]; |
||
| 1405 | ul[8] = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8]; |
||
| 1406 | |||
| 1407 | static double Wu[3]; |
||
| 1408 | if (nodeIOffset) { |
||
| 1409 | Wu[0] = nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5]; |
||
| 1410 | Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5]; |
||
| 1411 | Wu[2] = nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4]; |
||
| 1412 | |||
| 1413 | ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2]; |
||
| 1414 | ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 1415 | ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 1416 | } |
||
| 1417 | |||
| 1418 | if (nodeJOffset) { |
||
| 1419 | Wu[0] = nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11]; |
||
| 1420 | Wu[1] = -nodeJOffset[2]*ug[9] + nodeJOffset[0]*ug[11]; |
||
| 1421 | Wu[2] = nodeJOffset[1]*ug[9] - nodeJOffset[0]*ug[10]; |
||
| 1422 | |||
| 1423 | ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2]; |
||
| 1424 | ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2]; |
||
| 1425 | } |
||
| 1426 | |||
| 1427 | // compute displacements at point xi, in local coordinates |
||
| 1428 | static double uxl[3]; |
||
| 1429 | static Vector uxg(3); |
||
| 1430 | |||
| 1431 | uxl[0] = uxb(0) + ul[0]; |
||
| 1432 | uxl[1] = uxb(1) + (1-xi)*ul[1] + xi*ul[7]; |
||
| 1433 | uxl[2] = uxb(2) + (1-xi)*ul[2] + xi*ul[8]; |
||
| 1434 | |||
| 1435 | // rotate displacements to global coordinates |
||
| 1436 | // uxg = Rlj'*uxl |
||
| 1437 | //uxg.addMatrixTransposeVector(0.0, Rlj, uxl, 1.0); |
||
| 1438 | uxg(0) = R[0][0]*uxl[0] + R[1][0]*uxl[1] + R[2][0]*uxl[2]; |
||
| 1439 | uxg(1) = R[0][1]*uxl[0] + R[1][1]*uxl[1] + R[2][1]*uxl[2]; |
||
| 1440 | uxg(2) = R[0][2]*uxl[0] + R[1][2]*uxl[1] + R[2][2]*uxl[2]; |
||
| 1441 | |||
| 1442 | return uxg; |
||
| 272 | mhscott | 1443 | } |
| 1444 | |||
| 1445 | |||
| 1446 | void |
||
| 1271 | fmk | 1447 | PDeltaCrdTransf3d::Print(OPS_Stream &s, int flag) |
| 272 | mhscott | 1448 | { |
| 2221 | fmk | 1449 | s << "\nCrdTransf: " << this->getTag() << " Type: PDeltaCrdTransf3d" << endln; |
| 1450 | if (nodeIOffset) |
||
| 1451 | s << "\tNode I offset: " << nodeIOffset[0] << " " << nodeIOffset[1] << " "<< nodeIOffset[2] << endln; |
||
| 1452 | if (nodeJOffset) |
||
| 1453 | s << "\tNode J offset: " << nodeJOffset[0] << " " << nodeJOffset[1] << " "<< nodeJOffset[2] << endln; |
||
| 1454 | } |