My title

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.1 $
00022 // $Date: 2001/06/05 06:01:09 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/coordTransformation/PDeltaCrdTransf3d.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/crdTransf/PDeltaCrdTransf3d.C
00027 //
00028 // Written: Remo Magalhaes de Souza (rmsouza@ce.berkeley.edu)
00029 // Created: 04/2000
00030 // Revision: A
00031 // 
00032 // Purpose: This file contains the implementation for the 
00033 // PDeltaCrdTransf3d class. PDeltaCrdTransf3d 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 <iomanip.h>
00044 
00045 #include <PDeltaCrdTransf3d.h>
00046 
00047 // constructor:
00048 PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane):
00049   CrdTransf3d(tag, CRDTR_TAG_PDeltaCrdTransf3d),
00050   nodeIPtr(0), nodeJPtr(0),
00051   nodeIOffset(0), nodeJOffset(0),
00052   L(0), ul17(0), ul28(0)
00053 {
00054      R[2][0] = vecInLocXZPlane(0);
00055    R[2][1] = vecInLocXZPlane(1);
00056    R[2][2] = vecInLocXZPlane(2);
00057 
00058  // Does nothing
00059 }
00060 
00061 // constructor:
00062 PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane,
00063              const Vector &rigJntOffset1,
00064              const Vector &rigJntOffset2):
00065   CrdTransf3d(tag, CRDTR_TAG_PDeltaCrdTransf3d),
00066   nodeIPtr(0), nodeJPtr(0),
00067   nodeIOffset(0), nodeJOffset(0),
00068   L(0), ul17(0), ul28(0)
00069 {
00070    R[2][0] = vecInLocXZPlane(0);
00071    R[2][1] = vecInLocXZPlane(1);
00072    R[2][2] = vecInLocXZPlane(2);
00073 
00074  // check rigid joint offset for node I
00075  if (&rigJntOffset1 == 0 || rigJntOffset1.Size() != 3 ) {
00076   cerr << "PDeltaCrdTransf3d::PDeltaCrdTransf3d:  Invalid rigid joint offset vector for node I\n";
00077   cerr << "Size must be 3\n";      
00078  }
00079  else if (rigJntOffset1.Norm() > 0.0) {
00080   nodeIOffset = new double[3];
00081   nodeIOffset[0] = rigJntOffset1(0);
00082   nodeIOffset[1] = rigJntOffset1(1);
00083   nodeIOffset[2] = rigJntOffset1(2);
00084  }
00085    
00086    // check rigid joint offset for node J
00087  if (&rigJntOffset2 == 0 || rigJntOffset2.Size() != 3 ) {
00088   cerr << "PDeltaCrdTransf3d::PDeltaCrdTransf3d:  Invalid rigid joint offset vector for node J\n";
00089   cerr << "Size must be 3\n";      
00090  }
00091  else if (rigJntOffset2.Norm() > 0.0) {
00092   nodeJOffset = new double[3];
00093   nodeJOffset[0] = rigJntOffset2(0);
00094   nodeJOffset[1] = rigJntOffset2(1);
00095   nodeJOffset[2] = rigJntOffset2(2);
00096  }
00097 }
00098 
00099 
00100 
00101  
00102 // constructor:
00103 // invoked by a FEM_ObjectBroker, recvSelf() needs to be invoked on this object.
00104 PDeltaCrdTransf3d::PDeltaCrdTransf3d():
00105   CrdTransf3d(0, CRDTR_TAG_PDeltaCrdTransf3d),
00106   nodeIPtr(0), nodeJPtr(0),
00107   nodeIOffset(0), nodeJOffset(0),
00108   L(0), ul17(0), ul28(0)
00109 {
00110 
00111 }
00112 
00113 
00114 
00115 // destructor:
00116 PDeltaCrdTransf3d::~PDeltaCrdTransf3d() 
00117 {
00118  if (nodeIOffset)
00119   delete [] nodeIOffset;
00120  if (nodeJOffset)
00121   delete [] nodeJOffset;
00122 }
00123 
00124 
00125 int
00126 PDeltaCrdTransf3d::commitState(void)
00127 {
00128    return 0;
00129 }
00130 
00131 
00132 int
00133 PDeltaCrdTransf3d::revertToLastCommit(void)
00134 {
00135    return 0;
00136 }
00137 
00138 
00139 int
00140 PDeltaCrdTransf3d::revertToStart(void)
00141 {
00142    return 0;
00143 }
00144 
00145 
00146 int 
00147 PDeltaCrdTransf3d::initialize(Node *nodeIPointer, Node *nodeJPointer)
00148 {       
00149    int error;
00150 
00151    nodeIPtr = nodeIPointer;
00152    nodeJPtr = nodeJPointer;
00153 
00154    if ((!nodeIPtr) || (!nodeJPtr))
00155    {
00156       cerr << "\nPDeltaCrdTransf3d::initialize";
00157       cerr << "\ninvalid pointers to the element nodes\n";
00158       return -1;
00159    }
00160        
00161    // get element length and orientation
00162    if ((error = this->computeElemtLengthAndOrient()))
00163       return error;
00164       
00165    // get 3by3 rotation matrix
00166    if ((error = this->getLocalAxes()))
00167       return error;
00168 
00169    return 0;
00170 }
00171 
00172 
00173 int
00174 PDeltaCrdTransf3d::update(void)
00175 {
00176  const Vector &disp1 = nodeIPtr->getTrialDisp();
00177  const Vector &disp2 = nodeJPtr->getTrialDisp();
00178 
00179  static double ug[12];
00180  for (int i = 0; i < 6; i++) {
00181   ug[i]   = disp1(i);
00182   ug[i+6] = disp2(i);
00183  }
00184 
00185  double ul1, ul7, ul2, ul8;
00186 
00187  ul1 = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
00188  ul2 = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
00189 
00190  ul7 = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
00191  ul8 = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
00192 
00193  static double Wu[3];
00194 
00195  if (nodeIOffset) {
00196   Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
00197   Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
00198   Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
00199 
00200   ul1 += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00201   ul2 += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00202  }
00203 
00204  if (nodeJOffset) {
00205   Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
00206   Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
00207   Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
00208 
00209   ul7 += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00210   ul8 += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00211  }
00212 
00213  ul17 = ul1-ul7;
00214  ul28 = ul2-ul8;
00215 
00216  return 0;
00217 }
00218 
00219 
00220 int 
00221 PDeltaCrdTransf3d::computeElemtLengthAndOrient()
00222 {
00223    // element projection
00224    static Vector dx(3);
00225 
00226    const Vector &ndICoords = nodeIPtr->getCrds();
00227    const Vector &ndJCoords = nodeJPtr->getCrds();
00228 
00229    if (nodeIOffset == 0) {
00230     dx(0) = ndJCoords(0) - ndICoords(0);
00231     dx(1) = ndJCoords(1) - ndICoords(1);
00232     dx(2) = ndJCoords(2) - ndICoords(2);
00233    }
00234    else {
00235     dx(0) = ndJCoords(0) + nodeJOffset[0] - ndICoords(0) - nodeIOffset[0];
00236     dx(1) = ndJCoords(1) + nodeJOffset[1] - ndICoords(1) - nodeIOffset[1];
00237     dx(2) = ndJCoords(2) + nodeJOffset[2] - ndICoords(2) - nodeIOffset[2];
00238    }
00239    
00240    // calculate the element length
00241    L = dx.Norm();
00242 
00243    if (L == 0.0) {
00244       cerr << "\nPDeltaCrdTransf3d::computeElemtLengthAndOrien: 0 length\n";
00245       return -2;  
00246    }
00247 
00248    // calculate the element local x axis components (direction cossines)
00249    // wrt to the global coordinates
00250    R[0][0] = dx(0)/L;
00251    R[0][1] = dx(1)/L;
00252    R[0][2] = dx(2)/L;
00253    
00254    return 0;
00255 }
00256 
00257 
00258 int
00259 PDeltaCrdTransf3d::getLocalAxes(void)
00260 {
00261  // Compute y = v cross x
00262  // Note: v(i) is stored in R[2][i]
00263  static Vector vAxis(3);
00264  vAxis(0) = R[2][0]; vAxis(1) = R[2][1]; vAxis(2) = R[2][2];
00265 
00266  static Vector xAxis(3);
00267  xAxis(0) = R[0][0]; xAxis(1) = R[0][1]; xAxis(2) = R[0][2];
00268 
00269  static Vector yAxis(3);
00270 
00271  yAxis(0) = vAxis(1)*xAxis(2) - vAxis(2)*xAxis(1);
00272  yAxis(1) = vAxis(2)*xAxis(0) - vAxis(0)*xAxis(2);
00273  yAxis(2) = vAxis(0)*xAxis(1) - vAxis(1)*xAxis(0);
00274 
00275  double ynorm = yAxis.Norm();
00276 
00277  if (ynorm == 0) {
00278   cerr << "\nPDeltaCrdTransf3d::getLocalAxes";
00279   cerr << "\nvector v that defines plane xz is parallel to x axis\n";
00280   return -3;
00281  }
00282 
00283  yAxis /= ynorm;
00284 
00285  // Compute z = x cross y
00286  static Vector zAxis(3);
00287 
00288  zAxis(0) = xAxis(1)*yAxis(2) - xAxis(2)*yAxis(1);
00289  zAxis(1) = xAxis(2)*yAxis(0) - xAxis(0)*yAxis(2);
00290  zAxis(2) = xAxis(0)*yAxis(1) - xAxis(1)*yAxis(0);
00291 
00292  // Fill in transformation matrix
00293  R[1][0] = yAxis(0);
00294  R[1][1] = yAxis(1);
00295  R[1][2] = yAxis(2);
00296 
00297  R[2][0] = zAxis(0);
00298  R[2][1] = zAxis(1);
00299  R[2][2] = zAxis(2);
00300 
00301  return 0;
00302 }
00303 
00304 
00305 
00306 double 
00307 PDeltaCrdTransf3d::getInitialLength(void)
00308 {
00309    return L;
00310 }
00311 
00312 
00313 double 
00314 PDeltaCrdTransf3d::getDeformedLength(void)
00315 {
00316    return L;
00317 }
00318 
00319 
00320 const Vector &
00321 PDeltaCrdTransf3d::getBasicTrialDisp (void)
00322 {
00323  // determine global displacements
00324  const Vector &disp1 = nodeIPtr->getTrialDisp();
00325  const Vector &disp2 = nodeJPtr->getTrialDisp();
00326 
00327  static double ug[12];
00328  for (int i = 0; i < 6; i++) {
00329   ug[i]   = disp1(i);
00330   ug[i+6] = disp2(i);
00331  }
00332 
00333  double oneOverL = 1.0/L;
00334 
00335  static Vector ub(6);
00336 
00337  static double ul[12];
00338 
00339  ul[0]  = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2];
00340  ul[1]  = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
00341  ul[2]  = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
00342 
00343  ul[3]  = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5];
00344  ul[4]  = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5];
00345  ul[5]  = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5];
00346 
00347  ul[6]  = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8];
00348  ul[7]  = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
00349  ul[8]  = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
00350 
00351  ul[9]  = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11];
00352  ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11];
00353  ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11];
00354 
00355  static double Wu[3];
00356  if (nodeIOffset) {
00357   Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
00358   Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
00359   Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
00360 
00361   ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
00362   ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00363   ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00364  }
00365 
00366  if (nodeJOffset) {
00367   Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
00368   Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
00369   Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
00370 
00371   ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
00372   ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00373   ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00374  }
00375 
00376  ub(0) = ul[6] - ul[0];
00377  double tmp;
00378  tmp = oneOverL*(ul[1]-ul[7]);
00379  ub(1) = ul[5] + tmp;
00380  ub(2) = ul[11] + tmp;
00381  tmp = oneOverL*(ul[8]-ul[2]);
00382  ub(3) = ul[4] + tmp;
00383  ub(4) = ul[10] + tmp;
00384  ub(5) = ul[9] - ul[3];
00385 
00386  return ub;
00387 }
00388 
00389 
00390 const Vector &
00391 PDeltaCrdTransf3d::getBasicIncrDisp (void)
00392 {
00393  // determine global displacements
00394  const Vector &disp1 = nodeIPtr->getIncrDisp();
00395  const Vector &disp2 = nodeJPtr->getIncrDisp();
00396 
00397  static double ug[12];
00398  for (int i = 0; i < 6; i++) {
00399   ug[i]   = disp1(i);
00400   ug[i+6] = disp2(i);
00401  }
00402 
00403  double oneOverL = 1.0/L;
00404 
00405  static Vector ub(6);
00406 
00407  static double ul[12];
00408 
00409  ul[0]  = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2];
00410  ul[1]  = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
00411  ul[2]  = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
00412 
00413  ul[3]  = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5];
00414  ul[4]  = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5];
00415  ul[5]  = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5];
00416 
00417  ul[6]  = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8];
00418  ul[7]  = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
00419  ul[8]  = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
00420 
00421  ul[9]  = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11];
00422  ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11];
00423  ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11];
00424 
00425  static double Wu[3];
00426  if (nodeIOffset) {
00427   Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
00428   Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
00429   Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
00430 
00431   ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
00432   ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00433   ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00434  }
00435 
00436  if (nodeJOffset) {
00437   Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
00438   Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
00439   Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
00440 
00441   ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
00442   ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00443   ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00444  }
00445 
00446  ub(0) = ul[6] - ul[0];
00447  double tmp;
00448  tmp = oneOverL*(ul[1]-ul[7]);
00449  ub(1) = ul[5] + tmp;
00450  ub(2) = ul[11] + tmp;
00451  tmp = oneOverL*(ul[8]-ul[2]);
00452  ub(3) = ul[4] + tmp;
00453  ub(4) = ul[10] + tmp;
00454  ub(5) = ul[9] - ul[3];
00455 
00456  return ub;
00457 }
00458 
00459 
00460 const Vector &
00461 PDeltaCrdTransf3d::getBasicIncrDeltaDisp(void)
00462 {
00463  // determine global displacements
00464  const Vector &disp1 = nodeIPtr->getIncrDeltaDisp();
00465  const Vector &disp2 = nodeJPtr->getIncrDeltaDisp();
00466 
00467  static double ug[12];
00468  for (int i = 0; i < 6; i++) {
00469   ug[i]   = disp1(i);
00470   ug[i+6] = disp2(i);
00471  }
00472 
00473  double oneOverL = 1.0/L;
00474 
00475  static Vector ub(6);
00476 
00477  static double ul[12];
00478 
00479  ul[0]  = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2];
00480  ul[1]  = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
00481  ul[2]  = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
00482 
00483  ul[3]  = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5];
00484  ul[4]  = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5];
00485  ul[5]  = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5];
00486 
00487  ul[6]  = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8];
00488  ul[7]  = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
00489  ul[8]  = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
00490 
00491  ul[9]  = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11];
00492  ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11];
00493  ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11];
00494 
00495  static double Wu[3];
00496  if (nodeIOffset) {
00497   Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
00498   Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
00499   Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
00500 
00501   ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
00502   ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00503   ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00504  }
00505 
00506  if (nodeJOffset) {
00507   Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
00508   Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
00509   Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
00510 
00511   ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
00512   ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00513   ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00514  }
00515 
00516  ub(0) = ul[6] - ul[0];
00517  double tmp;
00518  tmp = oneOverL*(ul[1]-ul[7]);
00519  ub(1) = ul[5] + tmp;
00520  ub(2) = ul[11] + tmp;
00521  tmp = oneOverL*(ul[8]-ul[2]);
00522  ub(3) = ul[4] + tmp;
00523  ub(4) = ul[10] + tmp;
00524  ub(5) = ul[9] - ul[3];
00525 
00526  return ub;
00527 }
00528 
00529 
00530 const Vector &
00531 PDeltaCrdTransf3d::getGlobalResistingForce(const Vector &pb, const Vector &unifLoad)
00532 {
00533  // transform resisting forces from the basic system to local coordinates
00534  static double pl[12];
00535 
00536  double q0 = pb(0);
00537  double q1 = pb(1);
00538  double q2 = pb(2);
00539  double q3 = pb(3);
00540  double q4 = pb(4);
00541  double q5 = pb(5);
00542 
00543  double oneOverL = 1.0/L;
00544 
00545  pl[0]  = -q0;
00546  pl[1]  =  oneOverL*(q1+q2);
00547  pl[2]  = -oneOverL*(q3+q4);
00548  pl[3]  = -q5;
00549  pl[4]  =  q3;
00550  pl[5]  =  q1;
00551  pl[6]  =  q0;
00552  pl[7]  = -pl[1];
00553  pl[8]  = -pl[2];
00554  pl[9]  =  q5;
00555  pl[10] =  q4;
00556  pl[11] =  q2;
00557 
00558  // add end forces due to uniform distributed
00559  // loads to the system with rigid body modes
00560  pl[0] -= unifLoad(0)*L;
00561  double V;
00562  V = 0.5*unifLoad(1)*L;
00563  pl[1] -= V;
00564  pl[7] -= V;
00565  V = 0.5*unifLoad(2)*L;
00566  pl[2] -= V;
00567  pl[8] -= V;
00568 
00569  // Include leaning column effects (P-Delta)
00570  double NoverL;
00571  NoverL = ul17*q0*oneOverL;             
00572  pl[1] += NoverL;
00573  pl[7] -= NoverL;
00574  NoverL = ul28*q0*oneOverL;
00575  pl[2] += NoverL;
00576  pl[8] -= NoverL;
00577 
00578  // transform resisting forces  from local to global coordinates
00579  static Vector pg(12);
00580 
00581  pg(0)  = R[0][0]*pl[0] + R[1][0]*pl[1] + R[2][0]*pl[2];
00582  pg(1)  = R[0][1]*pl[0] + R[1][1]*pl[1] + R[2][1]*pl[2];
00583  pg(2)  = R[0][2]*pl[0] + R[1][2]*pl[1] + R[2][2]*pl[2];
00584 
00585  pg(3)  = R[0][0]*pl[3] + R[1][0]*pl[4] + R[2][0]*pl[5];
00586  pg(4)  = R[0][1]*pl[3] + R[1][1]*pl[4] + R[2][1]*pl[5];
00587  pg(5)  = R[0][2]*pl[3] + R[1][2]*pl[4] + R[2][2]*pl[5];
00588 
00589  pg(6)  = R[0][0]*pl[6] + R[1][0]*pl[7] + R[2][0]*pl[8];
00590  pg(7)  = R[0][1]*pl[6] + R[1][1]*pl[7] + R[2][1]*pl[8];
00591  pg(8)  = R[0][2]*pl[6] + R[1][2]*pl[7] + R[2][2]*pl[8];
00592 
00593  pg(9)  = R[0][0]*pl[9] + R[1][0]*pl[10] + R[2][0]*pl[11];
00594  pg(10) = R[0][1]*pl[9] + R[1][1]*pl[10] + R[2][1]*pl[11];
00595  pg(11) = R[0][2]*pl[9] + R[1][2]*pl[10] + R[2][2]*pl[11];
00596 
00597  if (nodeIOffset) {
00598   pg(3) += -nodeIOffset[2]*pg(1) + nodeIOffset[1]*pg(2);
00599   pg(4) +=  nodeIOffset[2]*pg(0) - nodeIOffset[0]*pg(2);
00600   pg(5) += -nodeIOffset[1]*pg(0) + nodeIOffset[0]*pg(1);
00601  }
00602 
00603  if (nodeJOffset) {
00604   pg(9)  += -nodeJOffset[2]*pg(7) + nodeJOffset[1]*pg(8);
00605   pg(10) +=  nodeJOffset[2]*pg(6) - nodeJOffset[0]*pg(8);
00606   pg(11) += -nodeJOffset[1]*pg(6) + nodeJOffset[0]*pg(7);
00607  }
00608  
00609  return pg;
00610 }
00611 
00612 const Matrix &
00613 PDeltaCrdTransf3d::getGlobalStiffMatrix (const Matrix &KB, const Vector &pb)
00614 {
00615  static Matrix kg(12,12); // Global stiffness for return
00616  static double kb[6][6];  // Basic stiffness
00617  static double kl[12][12]; // Local stiffness
00618  static double tmp[12][12]; // Temporary storage
00619 
00620  double oneOverL = 1.0/L;
00621 
00622  int i,j;
00623  for (i = 0; i < 6; i++)
00624   for (j = 0; j < 6; j++)
00625    kb[i][j] = KB(i,j);
00626 
00627  // Transform basic stiffness to local system
00628  // First compute kb*T_{bl}
00629  for (i = 0; i < 6; i++) {
00630   tmp[i][0]  = -kb[i][0];
00631   tmp[i][1]  =  oneOverL*(kb[i][1]+kb[i][2]);
00632   tmp[i][2]  = -oneOverL*(kb[i][3]+kb[i][4]);
00633   tmp[i][3]  = -kb[i][5];
00634   tmp[i][4]  =  kb[i][3];
00635   tmp[i][5]  =  kb[i][1];
00636   tmp[i][6]  =  kb[i][0];
00637   tmp[i][7]  = -tmp[i][1];
00638   tmp[i][8]  = -tmp[i][2];
00639   tmp[i][9]  =  kb[i][5];
00640   tmp[i][10] =  kb[i][4];
00641   tmp[i][11] =  kb[i][2];
00642  }
00643  
00644  // Now compute T'_{bl}*(kb*T_{bl})
00645  for (i = 0; i < 12; i++) {
00646   kl[0][i]  = -tmp[0][i];
00647   kl[1][i]  =  oneOverL*(tmp[1][i]+tmp[2][i]);
00648   kl[2][i]  = -oneOverL*(tmp[3][i]+tmp[4][i]);
00649   kl[3][i]  = -tmp[5][i];
00650   kl[4][i]  =  tmp[3][i];
00651   kl[5][i]  =  tmp[1][i];
00652   kl[6][i]  =  tmp[0][i];
00653   kl[7][i]  = -kl[1][i];
00654   kl[8][i]  = -kl[2][i];
00655   kl[9][i]  =  tmp[5][i];
00656   kl[10][i] =  tmp[4][i];
00657   kl[11][i] =  tmp[2][i];
00658  }
00659 
00660  // Include geometric stiffness effects in local system
00661  double NoverL = pb(0)*oneOverL;
00662  kl[1][1] += NoverL;
00663  kl[2][2] += NoverL;
00664  kl[7][7] += NoverL;
00665  kl[8][8] += NoverL;
00666  kl[1][7] -= NoverL;
00667  kl[7][1] -= NoverL;
00668  kl[2][8] -= NoverL;
00669  kl[8][2] -= NoverL;
00670  
00671  static double RWI[3][3];
00672 
00673  if (nodeIOffset) {
00674   // Compute RWI
00675   RWI[0][0] = -R[0][1]*nodeIOffset[2] + R[0][2]*nodeIOffset[1];
00676   RWI[1][0] = -R[1][1]*nodeIOffset[2] + R[1][2]*nodeIOffset[1];
00677   RWI[2][0] = -R[2][1]*nodeIOffset[2] + R[2][2]*nodeIOffset[1];
00678 
00679   RWI[0][1] =  R[0][0]*nodeIOffset[2] - R[0][2]*nodeIOffset[0];
00680   RWI[1][1] =  R[1][0]*nodeIOffset[2] - R[1][2]*nodeIOffset[0];
00681   RWI[2][1] =  R[2][0]*nodeIOffset[2] - R[2][2]*nodeIOffset[0];
00682 
00683   RWI[0][2] = -R[0][0]*nodeIOffset[1] + R[0][1]*nodeIOffset[0];
00684   RWI[1][2] = -R[1][0]*nodeIOffset[1] + R[1][1]*nodeIOffset[0];
00685   RWI[2][2] = -R[2][0]*nodeIOffset[1] + R[2][1]*nodeIOffset[0];
00686  }
00687 
00688  static double RWJ[3][3];
00689 
00690  if (nodeJOffset) {
00691   // Compute RWJ
00692   RWJ[0][0] = -R[0][1]*nodeJOffset[2] + R[0][2]*nodeJOffset[1];
00693   RWJ[1][0] = -R[1][1]*nodeJOffset[2] + R[1][2]*nodeJOffset[1];
00694   RWJ[2][0] = -R[2][1]*nodeJOffset[2] + R[2][2]*nodeJOffset[1];
00695 
00696   RWJ[0][1] =  R[0][0]*nodeJOffset[2] - R[0][2]*nodeJOffset[0];
00697   RWJ[1][1] =  R[1][0]*nodeJOffset[2] - R[1][2]*nodeJOffset[0];
00698   RWJ[2][1] =  R[2][0]*nodeJOffset[2] - R[2][2]*nodeJOffset[0];
00699 
00700   RWJ[0][2] = -R[0][0]*nodeJOffset[1] + R[0][1]*nodeJOffset[0];
00701   RWJ[1][2] = -R[1][0]*nodeJOffset[1] + R[1][1]*nodeJOffset[0];
00702   RWJ[2][2] = -R[2][0]*nodeJOffset[1] + R[2][1]*nodeJOffset[0];
00703  }
00704 
00705  // Transform local stiffness to global system
00706  // First compute kl*T_{lg}
00707  int m;
00708  for (m = 0; m < 12; m++) {
00709   tmp[m][0] = kl[m][0]*R[0][0] + kl[m][1]*R[1][0]  + kl[m][2]*R[2][0];
00710   tmp[m][1] = kl[m][0]*R[0][1] + kl[m][1]*R[1][1]  + kl[m][2]*R[2][1];
00711   tmp[m][2] = kl[m][0]*R[0][2] + kl[m][1]*R[1][2]  + kl[m][2]*R[2][2];
00712 
00713   tmp[m][3] = kl[m][3]*R[0][0] + kl[m][4]*R[1][0]  + kl[m][5]*R[2][0];
00714   tmp[m][4] = kl[m][3]*R[0][1] + kl[m][4]*R[1][1]  + kl[m][5]*R[2][1];
00715   tmp[m][5] = kl[m][3]*R[0][2] + kl[m][4]*R[1][2]  + kl[m][5]*R[2][2];
00716 
00717   if (nodeIOffset) {
00718    tmp[m][3]  += kl[m][0]*RWI[0][0]  + kl[m][1]*RWI[1][0]  + kl[m][2]*RWI[2][0];
00719    tmp[m][4]  += kl[m][0]*RWI[0][1]  + kl[m][1]*RWI[1][1]  + kl[m][2]*RWI[2][1];
00720    tmp[m][5]  += kl[m][0]*RWI[0][2]  + kl[m][1]*RWI[1][2]  + kl[m][2]*RWI[2][2];
00721   }
00722 
00723   tmp[m][6] = kl[m][6]*R[0][0] + kl[m][7]*R[1][0]  + kl[m][8]*R[2][0];
00724   tmp[m][7] = kl[m][6]*R[0][1] + kl[m][7]*R[1][1]  + kl[m][8]*R[2][1];
00725   tmp[m][8] = kl[m][6]*R[0][2] + kl[m][7]*R[1][2]  + kl[m][8]*R[2][2];
00726 
00727   tmp[m][9]  = kl[m][9]*R[0][0] + kl[m][10]*R[1][0] + kl[m][11]*R[2][0];
00728   tmp[m][10] = kl[m][9]*R[0][1] + kl[m][10]*R[1][1] + kl[m][11]*R[2][1];
00729   tmp[m][11] = kl[m][9]*R[0][2] + kl[m][10]*R[1][2] + kl[m][11]*R[2][2];
00730 
00731   if (nodeJOffset) {
00732    tmp[m][9]   += kl[m][6]*RWJ[0][0]  + kl[m][7]*RWJ[1][0]  + kl[m][8]*RWJ[2][0];
00733    tmp[m][10]  += kl[m][6]*RWJ[0][1]  + kl[m][7]*RWJ[1][1]  + kl[m][8]*RWJ[2][1];
00734    tmp[m][11]  += kl[m][6]*RWJ[0][2]  + kl[m][7]*RWJ[1][2]  + kl[m][8]*RWJ[2][2];
00735   }
00736 
00737  }
00738 
00739  // Now compute T'_{lg}*(kl*T_{lg})
00740  for (m = 0; m < 12; m++) {
00741   kg(0,m) = R[0][0]*tmp[0][m] + R[1][0]*tmp[1][m]  + R[2][0]*tmp[2][m];
00742   kg(1,m) = R[0][1]*tmp[0][m] + R[1][1]*tmp[1][m]  + R[2][1]*tmp[2][m];
00743   kg(2,m) = R[0][2]*tmp[0][m] + R[1][2]*tmp[1][m]  + R[2][2]*tmp[2][m];
00744 
00745   kg(3,m) = R[0][0]*tmp[3][m] + R[1][0]*tmp[4][m]  + R[2][0]*tmp[5][m];
00746   kg(4,m) = R[0][1]*tmp[3][m] + R[1][1]*tmp[4][m]  + R[2][1]*tmp[5][m];
00747   kg(5,m) = R[0][2]*tmp[3][m] + R[1][2]*tmp[4][m]  + R[2][2]*tmp[5][m];
00748 
00749   if (nodeIOffset) {
00750    kg(3,m) += RWI[0][0]*tmp[0][m]  + RWI[1][0]*tmp[1][m] + RWI[2][0]*tmp[2][m];
00751    kg(4,m) += RWI[0][1]*tmp[0][m]  + RWI[1][1]*tmp[1][m] + RWI[2][1]*tmp[2][m];
00752    kg(5,m) += RWI[0][2]*tmp[0][m]  + RWI[1][2]*tmp[1][m] + RWI[2][2]*tmp[2][m];
00753   }
00754 
00755   kg(6,m) = R[0][0]*tmp[6][m] + R[1][0]*tmp[7][m]  + R[2][0]*tmp[8][m];
00756   kg(7,m) = R[0][1]*tmp[6][m] + R[1][1]*tmp[7][m]  + R[2][1]*tmp[8][m];
00757   kg(8,m) = R[0][2]*tmp[6][m] + R[1][2]*tmp[7][m]  + R[2][2]*tmp[8][m];
00758 
00759   kg(9,m)  = R[0][0]*tmp[9][m] + R[1][0]*tmp[10][m] + R[2][0]*tmp[11][m];
00760   kg(10,m) = R[0][1]*tmp[9][m] + R[1][1]*tmp[10][m] + R[2][1]*tmp[11][m];
00761   kg(11,m) = R[0][2]*tmp[9][m] + R[1][2]*tmp[10][m] + R[2][2]*tmp[11][m];
00762 
00763   if (nodeJOffset) {
00764    kg(9,m)  += RWJ[0][0]*tmp[6][m]  + RWJ[1][0]*tmp[7][m] + RWJ[2][0]*tmp[8][m];
00765    kg(10,m) += RWJ[0][1]*tmp[6][m]  + RWJ[1][1]*tmp[7][m] + RWJ[2][1]*tmp[8][m];
00766    kg(11,m) += RWJ[0][2]*tmp[6][m]  + RWJ[1][2]*tmp[7][m] + RWJ[2][2]*tmp[8][m];
00767   }
00768  }
00769 
00770  return kg;
00771 }
00772   
00773 
00774 
00775 
00776 CrdTransf3d *
00777 PDeltaCrdTransf3d::getCopy(void)
00778 {
00779   // create a new instance of PDeltaCrdTransf3d 
00780 
00781   PDeltaCrdTransf3d *theCopy;
00782 
00783   static Vector xz(3);
00784   xz(0) = R[2][0];
00785   xz(1) = R[2][1];
00786   xz(2) = R[2][2];
00787 
00788   Vector offsetI(3);
00789   Vector offsetJ(3);
00790 
00791   if (nodeIOffset) {
00792    offsetI(0) = nodeIOffset[0];
00793    offsetI(1) = nodeIOffset[1];
00794    offsetI(2) = nodeIOffset[2];
00795   }
00796 
00797   if (nodeJOffset) {
00798    offsetJ(0) = nodeJOffset[0];
00799    offsetJ(1) = nodeJOffset[1];
00800    offsetJ(2) = nodeJOffset[2];
00801   }
00802   
00803   theCopy = new PDeltaCrdTransf3d(this->getTag(), xz, offsetI, offsetJ);
00804 
00805   theCopy->nodeIPtr = nodeIPtr;
00806   theCopy->nodeJPtr = nodeJPtr;
00807   theCopy->L = L;
00808   theCopy->ul17 = ul17;
00809   theCopy->ul28 = ul28;
00810   for (int i = 0; i < 3; i++)
00811    for (int j = 0; j < 3; j++)
00812     theCopy->R[i][j] = R[i][j];
00813 
00814   
00815   return theCopy;
00816 }
00817 
00818 
00819 int 
00820 PDeltaCrdTransf3d::sendSelf(int cTag, Channel &theChannel)
00821 {
00822  int res = 0;
00823 
00824  static Vector data(9);
00825 
00826  data(0) = this->getTag();
00827  data(6) = L;
00828 
00829  res += theChannel.sendVector(this->getDbTag(), cTag, data);
00830  if (res < 0) {
00831   g3ErrorHandler->warning("%s - failed to send Vector",
00832    "PDeltaCrdTransf3d::sendSelf");
00833   return res;
00834  }
00835 
00836     return res;
00837 }
00838 
00839     
00840 
00841 int 
00842 PDeltaCrdTransf3d::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
00843 {
00844  int res = 0;
00845 
00846  static Vector data(9);
00847 
00848  res += theChannel.recvVector(this->getDbTag(), cTag, data);
00849  if (res < 0) {
00850   g3ErrorHandler->warning("%s - failed to receive Vector",
00851    "PDeltaCrdTransf3d::recvSelf");
00852   return res;
00853  }
00854 
00855  this->setTag((int)data(0));
00856  L = data(6);
00857 
00858     return res;
00859 }
00860   
00861 
00862 const Vector &
00863 PDeltaCrdTransf3d::getPointGlobalCoordFromLocal(const Vector &xl)
00864 {
00865    static Vector xg(3);
00866 
00867    //xg = nodeIPtr->getCrds() + nodeIOffset;
00868    xg = nodeIPtr->getCrds();
00869    
00870    if (nodeIOffset) {
00871     xg(0) += nodeIOffset[0];
00872     xg(1) += nodeIOffset[1];
00873     xg(2) += nodeIOffset[2];
00874    }
00875 
00876    // xg = xg + Rlj'*xl
00877    //xg.addMatrixTransposeVector(1.0, Rlj, xl, 1.0);
00878    xg(0) += R[0][0]*xl(0) + R[0][1]*xl(1) + R[0][2]*xl(2);
00879    xg(1) += R[1][0]*xl(0) + R[1][1]*xl(1) + R[1][2]*xl(2);
00880    xg(2) += R[2][0]*xl(0) + R[2][1]*xl(1) + R[2][2]*xl(2);
00881      
00882    return xg;  
00883 }
00884 
00885     
00886 const Vector &
00887 PDeltaCrdTransf3d::getPointGlobalDisplFromBasic (double xi, const Vector &uxb)
00888 {
00889    // determine global displacements
00890    const Vector &disp1 = nodeIPtr->getTrialDisp();
00891    const Vector &disp2 = nodeJPtr->getTrialDisp();
00892 
00893    static double ug[12];
00894    for (int i = 0; i < 6; i++)
00895    {
00896       ug[i]   = disp1(i);
00897       ug[i+6] = disp2(i);
00898    }
00899 
00900    // transform global end displacements to local coordinates
00901    //ul.addMatrixVector(0.0, Tlg,  ug, 1.0);       //  ul = Tlg *  ug;
00902  static double ul[12];
00903 
00904  ul[0]  = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2];
00905  ul[1]  = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
00906  ul[2]  = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
00907 
00908  ul[7]  = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
00909  ul[8]  = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
00910 
00911  static double Wu[3];
00912  if (nodeIOffset) {
00913   Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
00914   Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
00915   Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
00916 
00917   ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
00918   ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00919   ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00920  }
00921 
00922  if (nodeJOffset) {
00923   Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
00924   Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
00925   Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
00926 
00927   ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
00928   ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
00929  }
00930  
00931    // compute displacements at point xi, in local coordinates
00932    static double uxl[3];
00933    static Vector uxg(3);
00934 
00935    uxl[0] = uxb(0) +        ul[0];
00936    uxl[1] = uxb(1) + (1-xi)*ul[1] + xi*ul[7];
00937    uxl[2] = uxb(2) + (1-xi)*ul[2] + xi*ul[8];
00938   
00939    // rotate displacements to global coordinates
00940    // uxg = Rlj'*uxl
00941    //uxg.addMatrixTransposeVector(0.0, Rlj, uxl, 1.0);
00942    uxg(0) = R[0][0]*uxl[0] + R[1][0]*uxl[1] + R[2][0]*uxl[2];
00943    uxg(1) = R[0][1]*uxl[0] + R[1][1]*uxl[1] + R[2][1]*uxl[2];
00944    uxg(2) = R[0][2]*uxl[0] + R[1][2]*uxl[1] + R[2][2]*uxl[2];
00945      
00946    return uxg;  
00947 }
00948 
00949 
00950 
00951 void
00952 PDeltaCrdTransf3d::Print(ostream &s, int flag)
00953 {
00954    s << "\nCrdTransf: " << this->getTag() << " Type: PDeltaCrdTransf3d" << endl;
00955    if (nodeIOffset)
00956     s << "\tNode I offset: " << nodeIOffset[0] << ' ' << nodeIOffset[1] << ' '<< nodeIOffset[2] << endl;
00957    if (nodeJOffset)
00958     s << "\tNode J offset: " << nodeJOffset[0] << ' ' << nodeJOffset[1] << ' '<< nodeJOffset[2] << endl;
00959 }
00960 
00961 
00962 
00963 
00964 
00965 
00966 
00967