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