Subversion Repositories OpenSees

Rev

Rev 1894 | Rev 2891 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 1894 Rev 2221
Line 15... Line 15...
15
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
15
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
16
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
16
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
17
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
17
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
18
**                                                                    **
18
**                                                                    **
19
** ****************************************************************** */
19
** ****************************************************************** */
20
                                                                       
-
 
21
// $Revision: 1.12 $
-
 
22
// $Date: 2004-10-08 22:03:04 $
-
 
-
 
20
-
 
21
// $Revision: 1.13 $
-
 
22
// $Date: 2005-12-15 00:30:38 $
23
// $Source: /usr/local/cvs/OpenSees/SRC/coordTransformation/PDeltaCrdTransf3d.cpp,v $
23
// $Source: /usr/local/cvs/OpenSees/SRC/coordTransformation/PDeltaCrdTransf3d.cpp,v $
24
                                                                       
-
 
25
                                                                       
-
 
-
 
24
-
 
25
26
// Written: Remo Magalhaes de Souza (rmsouza@ce.berkeley.edu)
26
// Written: Remo Magalhaes de Souza (rmsouza@ce.berkeley.edu)
27
// Created: 04/2000
27
// Created: 04/2000
28
// Revision: A
28
// Revision: A
-
 
29
//
-
 
30
// Modified: 04/2005 Andreas Schellenberg (getBasicTrialVel, getBasicTrialAccel)
29
// 
31
// 
30
// Purpose: This file contains the implementation for the 
32
// Purpose: This file contains the implementation for the 
31
// PDeltaCrdTransf3d class. PDeltaCrdTransf3d is a linear
33
// PDeltaCrdTransf3d class. PDeltaCrdTransf3d is a linear
32
// transformation for a planar frame between the global 
34
// transformation for a planar frame between the global 
33
// and basic coordinate systems
35
// and basic coordinate systems
Line 41... Line 43...
41
#include <PDeltaCrdTransf3d.h>
43
#include <PDeltaCrdTransf3d.h>
42
44
43
45
44
// constructor:
46
// constructor:
45
PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane):
47
PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane):
46
  CrdTransf3d(tag, CRDTR_TAG_PDeltaCrdTransf3d),
-
 
47
  nodeIPtr(0), nodeJPtr(0),
-
 
48
  nodeIOffset(0), nodeJOffset(0),
-
 
49
  L(0), ul17(0), ul28(0),
-
 
50
  nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false)
-
 
-
 
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)
51
{
53
{
52
        for (int i = 0; i < 2; i++)
-
 
53
                for (int j = 0; j < 3; j++)
-
 
54
                        R[i][j] = 0.0;
-
 
-
 
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
}
55
64
56
        R[2][0] = vecInLocXZPlane(0);
-
 
57
        R[2][1] = vecInLocXZPlane(1);
-
 
58
        R[2][2] = vecInLocXZPlane(2);
-
 
59
-
 
60
        // Does nothing
-
 
61
}
-
 
62
65
63
// constructor:
66
// constructor:
64
PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane,
67
PDeltaCrdTransf3d::PDeltaCrdTransf3d(int tag, const Vector &vecInLocXZPlane,
65
                                     const Vector &rigJntOffset1,
-
 
66
                                     const Vector &rigJntOffset2):
-
 
67
  CrdTransf3d(tag, CRDTR_TAG_PDeltaCrdTransf3d),
-
 
68
  nodeIPtr(0), nodeJPtr(0),
-
 
69
  nodeIOffset(0), nodeJOffset(0),
-
 
70
  L(0), ul17(0), ul28(0),
-
 
71
  nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false)
-
 
-
 
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)
72
{
75
{
73
        for (int i = 0; i < 2; i++)
-
 
74
                for (int j = 0; j < 3; j++)
-
 
75
                        R[i][j] = 0.0;
-
 
76
-
 
77
        R[2][0] = vecInLocXZPlane(0);
-
 
78
        R[2][1] = vecInLocXZPlane(1);
-
 
79
        R[2][2] = vecInLocXZPlane(2);
-
 
80
-
 
81
        // check rigid joint offset for node I
-
 
82
        if (&rigJntOffset1 == 0 || rigJntOffset1.Size() != 3 ) {
-
 
83
                opserr << "PDeltaCrdTransf3d::PDeltaCrdTransf3d:  Invalid rigid joint offset vector for node I\n";
-
 
84
                opserr << "Size must be 3\n";      
-
 
85
        }
-
 
86
        else if (rigJntOffset1.Norm() > 0.0) {
-
 
87
                nodeIOffset = new double[3];
-
 
88
                nodeIOffset[0] = rigJntOffset1(0);
-
 
89
                nodeIOffset[1] = rigJntOffset1(1);
-
 
90
                nodeIOffset[2] = rigJntOffset1(2);
-
 
91
        }
-
 
92
   
-
 
93
   // check rigid joint offset for node J
-
 
94
        if (&rigJntOffset2 == 0 || rigJntOffset2.Size() != 3 ) {
-
 
95
                opserr << "PDeltaCrdTransf3d::PDeltaCrdTransf3d:  Invalid rigid joint offset vector for node J\n";
-
 
96
                opserr << "Size must be 3\n";      
-
 
97
        }
-
 
98
        else if (rigJntOffset2.Norm() > 0.0) {
-
 
99
                nodeJOffset = new double[3];
-
 
100
                nodeJOffset[0] = rigJntOffset2(0);
-
 
101
                nodeJOffset[1] = rigJntOffset2(1);
-
 
102
                nodeJOffset[2] = rigJntOffset2(2);
-
 
103
        }
-
 
-
 
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
        }
104
}
107
}
105
108
106
109
107
-
 
108
 
-
 
109
// constructor:
110
// constructor:
110
// invoked by a FEM_ObjectBroker, recvSelf() needs to be invoked on this object.
111
// invoked by a FEM_ObjectBroker, recvSelf() needs to be invoked on this object.
111
PDeltaCrdTransf3d::PDeltaCrdTransf3d():
112
PDeltaCrdTransf3d::PDeltaCrdTransf3d():
112
  CrdTransf3d(0, CRDTR_TAG_PDeltaCrdTransf3d),
-
 
113
  nodeIPtr(0), nodeJPtr(0),
-
 
114
  nodeIOffset(0), nodeJOffset(0),
-
 
115
  L(0), ul17(0), ul28(0),
-
 
116
  nodeIInitialDisp(0), nodeJInitialDisp(0), initialDispChecked(false)
-
 
-
 
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)
117
{
118
{
118
        for (int i = 0; i < 3; i++)
-
 
119
                for (int j = 0; j < 3; j++)
-
 
120
                        R[i][j] = 0.0;
-
 
-
 
119
    for (int i = 0; i < 3; i++)
-
 
120
        for (int j = 0; j < 3; j++)
-
 
121
            R[i][j] = 0.0;
121
}
122
}
122
-
 
123
123
124
124
125
// destructor:
125
// destructor:
126
PDeltaCrdTransf3d::~PDeltaCrdTransf3d()
126
PDeltaCrdTransf3d::~PDeltaCrdTransf3d()
127
{
127
{
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;
-
 
-
 
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;
136
}
136
}
137
137
138
138
139
int
139
int
140
PDeltaCrdTransf3d::commitState(void)
140
PDeltaCrdTransf3d::commitState(void)
141
{
141
{
142
   return 0;
-
 
-
 
142
    return 0;
143
}
143
}
144
144
145
145
146
int
146
int
147
PDeltaCrdTransf3d::revertToLastCommit(void)
147
PDeltaCrdTransf3d::revertToLastCommit(void)
148
{
148
{
149
   return 0;
-
 
-
 
149
    return 0;
150
}
150
}
151
151
152
152
153
int
153
int
154
PDeltaCrdTransf3d::revertToStart(void)
154
PDeltaCrdTransf3d::revertToStart(void)
155
{
155
{
156
   return 0;
-
 
-
 
156
    return 0;
157
}
157
}
158
158
159
159
160
int
160
int
161
PDeltaCrdTransf3d::initialize(Node *nodeIPointer, Node *nodeJPointer)
161
PDeltaCrdTransf3d::initialize(Node *nodeIPointer, Node *nodeJPointer)
162
{      
162
{      
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;
-
 
-
 
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;
211
}
211
}
212
212
213
213
214
int
214
int
215
PDeltaCrdTransf3d::update(void)
215
PDeltaCrdTransf3d::update(void)
216
{
216
{
217
  const Vector &disp1 = nodeIPtr->getTrialDisp();
-
 
218
  const Vector &disp2 = nodeJPtr->getTrialDisp();
-
 
219
-
 
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
  }
-
 
225
-
 
226
  if (nodeIInitialDisp != 0) {
-
 
227
    for (int j=0; j<6; j++)
-
 
228
      ug[j] -= nodeIInitialDisp[j];
-
 
229
  }
-
 
-
 
217
    const Vector &disp1 = nodeIPtr->getTrialDisp();
-
 
218
    const Vector &disp2 = nodeJPtr->getTrialDisp();
-
 
219
   
-
 
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
    }
-
 
225
   
-
 
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];
230
   
243
   
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];
-
 
-
 
244
    static double Wu[3];
250
   
245
   
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;
-
 
-
 
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;
268
}
268
}
269
269
270
270
271
int
271
int
272
PDeltaCrdTransf3d::computeElemtLengthAndOrient()
272
PDeltaCrdTransf3d::computeElemtLengthAndOrient()
273
{
273
{
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;
-
 
-
 
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;
323
}
323
}
324
324
325
325
326
int
326
int
327
PDeltaCrdTransf3d::getLocalAxes(Vector &XAxis, Vector &YAxis, Vector &ZAxis)
327
PDeltaCrdTransf3d::getLocalAxes(Vector &XAxis, Vector &YAxis, Vector &ZAxis)
328
{
328
{
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;
-
 
-
 
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;
374
}
374
}
375
-
 
376
375
377
376
378
double
377
double
379
PDeltaCrdTransf3d::getInitialLength(void)
378
PDeltaCrdTransf3d::getInitialLength(void)
380
{
379
{
381
   return L;
-
 
-
 
380
    return L;
382
}
381
}
383
382
384
383
385
double
384
double
386
PDeltaCrdTransf3d::getDeformedLength(void)
385
PDeltaCrdTransf3d::getDeformedLength(void)
387
{
386
{
388
   return L;
-
 
-
 
387
    return L;
389
}
388
}
390
389
391
390
392
const Vector &
391
const Vector &
393
PDeltaCrdTransf3d::getBasicTrialDisp (void)
392
PDeltaCrdTransf3d::getBasicTrialDisp (void)
394
{
393
{
395
  // determine global displacements
-
 
396
  const Vector &disp1 = nodeIPtr->getTrialDisp();
-
 
397
  const Vector &disp2 = nodeJPtr->getTrialDisp();
-
 
398
 
-
 
399
  static double ug[12];
-
 
400
  for (int i = 0; i < 6; i++) {
-
 
401
    ug[i]   = disp1(i);
-
 
402
    ug[i+6] = disp2(i);
-
 
403
  }
-
 
404
 
-
 
405
  if (nodeIInitialDisp != 0) {
-
 
406
    for (int j=0; j<6; j++)
-
 
407
      ug[j] -= nodeIInitialDisp[j];
-
 
408
  }
-
 
409
 
-
 
410
  if (nodeJInitialDisp != 0) {
-
 
411
    for (int j=0; j<6; j++)
-
 
412
      ug[j+6] -= nodeJInitialDisp[j];
-
 
413
  }
-
 
414
 
-
 
415
  double oneOverL = 1.0/L;
-
 
416
 
-
 
417
  static Vector ub(6);
-
 
418
 
-
 
419
  static double ul[12];
-
 
420
 
-
 
421
  ul[0]  = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2];
-
 
422
  ul[1]  = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
-
 
423
  ul[2]  = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
-
 
424
 
-
 
425
  ul[3]  = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5];
-
 
426
  ul[4]  = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5];
-
 
427
  ul[5]  = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5];
-
 
428
 
-
 
429
  ul[6]  = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8];
-
 
430
  ul[7]  = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
-
 
431
  ul[8]  = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
-
 
432
 
-
 
433
  ul[9]  = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11];
-
 
434
  ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11];
-
 
435
  ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11];
-
 
436
 
-
 
437
  static double Wu[3];
-
 
438
  if (nodeIOffset) {
-
 
439
    Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
-
 
440
    Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
-
 
441
    Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
-
 
-
 
394
    // determine global displacements
-
 
395
    const Vector &disp1 = nodeIPtr->getTrialDisp();
-
 
396
    const Vector &disp2 = nodeJPtr->getTrialDisp();
442
   
397
   
443
    ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
-
 
444
    ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
-
 
445
    ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
-
 
446
  }
-
 
447
 
-
 
448
  if (nodeJOffset) {
-
 
449
    Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
-
 
450
    Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
-
 
451
    Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
-
 
-
 
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
    }
452
   
403
   
453
    ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
-
 
454
    ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
-
 
455
    ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
-
 
456
  }
-
 
457
 
-
 
458
  ub(0) = ul[6] - ul[0];
-
 
459
  double tmp;
-
 
460
  tmp = oneOverL*(ul[1]-ul[7]);
-
 
461
  ub(1) = ul[5] + tmp;
-
 
462
  ub(2) = ul[11] + tmp;
-
 
463
  tmp = oneOverL*(ul[8]-ul[2]);
-
 
464
  ub(3) = ul[4] + tmp;
-
 
465
  ub(4) = ul[10] + tmp;
-
 
466
  ub(5) = ul[9] - ul[3];
-
 
467
 
-
 
468
  return ub;
-
 
-
 
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;
469
}
468
}
470
469
471
470
472
const Vector &
471
const Vector &
473
PDeltaCrdTransf3d::getBasicIncrDisp (void)
472
PDeltaCrdTransf3d::getBasicIncrDisp (void)
474
{
473
{
475
        // determine global displacements
-
 
476
        const Vector &disp1 = nodeIPtr->getIncrDisp();
-
 
477
        const Vector &disp2 = nodeJPtr->getIncrDisp();
-
 
-
 
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
}
478
539
479
        static double ug[12];
-
 
480
        for (int i = 0; i < 6; i++) {
-
 
481
                ug[i]   = disp1(i);
-
 
482
                ug[i+6] = disp2(i);
-
 
483
        }
-
 
484
540
485
        double oneOverL = 1.0/L;
-
 
-
 
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
}
486
609
487
        static Vector ub(6);
-
 
488
-
 
489
        static double ul[12];
-
 
490
-
 
491
        ul[0]  = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2];
-
 
492
        ul[1]  = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
-
 
493
        ul[2]  = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
-
 
494
-
 
495
        ul[3]  = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5];
-
 
496
        ul[4]  = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5];
-
 
497
        ul[5]  = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5];
-
 
498
-
 
499
        ul[6]  = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8];
-
 
500
        ul[7]  = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
-
 
501
        ul[8]  = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
-
 
502
-
 
503
        ul[9]  = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11];
-
 
504
        ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11];
-
 
505
        ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11];
-
 
506
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];
-
 
619
        for (int i = 0; i < 6; i++) {
-
 
620
                vg[i]   = vel1(i);
-
 
621
                vg[i+6] = vel2(i);
-
 
622
        }
-
 
623
       
-
 
624
        double oneOverL = 1.0/L;
-
 
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
       
507
        static double Wu[3];
646
        static double Wu[3];
508
        if (nodeIOffset) {
647
        if (nodeIOffset) {
509
                Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
-
 
510
                Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
-
 
511
                Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
-
 
512
-
 
513
                ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
-
 
514
                ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
-
 
515
                ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
-
 
-
 
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];
516
        }
655
        }
517
-
 
-
 
656
       
518
        if (nodeJOffset) {
657
        if (nodeJOffset) {
519
                Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
-
 
520
                Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
-
 
521
                Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
-
 
522
-
 
523
                ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
-
 
524
                ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
-
 
525
                ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
-
 
-
 
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];
526
        }
665
        }
527
-
 
528
        ub(0) = ul[6] - ul[0];
-
 
-
 
666
       
-
 
667
        vb(0) = vl[6] - vl[0];
529
        double tmp;
668
        double tmp;
530
        tmp = oneOverL*(ul[1]-ul[7]);
-
 
531
        ub(1) = ul[5] + tmp;
-
 
532
        ub(2) = ul[11] + tmp;
-
 
533
        tmp = oneOverL*(ul[8]-ul[2]);
-
 
534
        ub(3) = ul[4] + tmp;
-
 
535
        ub(4) = ul[10] + tmp;
-
 
536
        ub(5) = ul[9] - ul[3];
-
 
537
-
 
538
        return ub;
-
 
-
 
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;
539
}
678
}
540
679
541
680
542
const Vector &
681
const Vector &
543
PDeltaCrdTransf3d::getBasicIncrDeltaDisp(void)
-
 
-
 
682
PDeltaCrdTransf3d::getBasicTrialAccel(void)
544
{
683
{
545
        // determine global displacements
-
 
546
        const Vector &disp1 = nodeIPtr->getIncrDeltaDisp();
-
 
547
        const Vector &disp2 = nodeJPtr->getIncrDeltaDisp();
-
 
548
-
 
549
        static double ug[12];
-
 
-
 
684
        // determine global accelerations
-
 
685
        const Vector &accel1 = nodeIPtr->getTrialAccel();
-
 
686
        const Vector &accel2 = nodeJPtr->getTrialAccel();
-
 
687
       
-
 
688
        static double ag[12];
550
        for (int i = 0; i < 6; i++) {
689
        for (int i = 0; i < 6; i++) {
551
                ug[i]   = disp1(i);
-
 
552
                ug[i+6] = disp2(i);
-
 
-
 
690
                ag[i]   = accel1(i);
-
 
691
                ag[i+6] = accel2(i);
553
        }
692
        }
554
-
 
-
 
693
       
555
        double oneOverL = 1.0/L;
694
        double oneOverL = 1.0/L;
556
-
 
557
        static Vector ub(6);
-
 
558
-
 
559
        static double ul[12];
-
 
560
-
 
561
        ul[0]  = R[0][0]*ug[0] + R[0][1]*ug[1] + R[0][2]*ug[2];
-
 
562
        ul[1]  = R[1][0]*ug[0] + R[1][1]*ug[1] + R[1][2]*ug[2];
-
 
563
        ul[2]  = R[2][0]*ug[0] + R[2][1]*ug[1] + R[2][2]*ug[2];
-
 
564
-
 
565
        ul[3]  = R[0][0]*ug[3] + R[0][1]*ug[4] + R[0][2]*ug[5];
-
 
566
        ul[4]  = R[1][0]*ug[3] + R[1][1]*ug[4] + R[1][2]*ug[5];
-
 
567
        ul[5]  = R[2][0]*ug[3] + R[2][1]*ug[4] + R[2][2]*ug[5];
-
 
568
-
 
569
        ul[6]  = R[0][0]*ug[6] + R[0][1]*ug[7] + R[0][2]*ug[8];
-
 
570
        ul[7]  = R[1][0]*ug[6] + R[1][1]*ug[7] + R[1][2]*ug[8];
-
 
571
        ul[8]  = R[2][0]*ug[6] + R[2][1]*ug[7] + R[2][2]*ug[8];
-
 
572
-
 
573
        ul[9]  = R[0][0]*ug[9] + R[0][1]*ug[10] + R[0][2]*ug[11];
-
 
574
        ul[10] = R[1][0]*ug[9] + R[1][1]*ug[10] + R[1][2]*ug[11];
-
 
575
        ul[11] = R[2][0]*ug[9] + R[2][1]*ug[10] + R[2][2]*ug[11];
-
 
576
-
 
-
 
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
       
577
        static double Wu[3];
716
        static double Wu[3];
578
        if (nodeIOffset) {
717
        if (nodeIOffset) {
579
                Wu[0] =  nodeIOffset[2]*ug[4] - nodeIOffset[1]*ug[5];
-
 
580
                Wu[1] = -nodeIOffset[2]*ug[3] + nodeIOffset[0]*ug[5];
-
 
581
                Wu[2] =  nodeIOffset[1]*ug[3] - nodeIOffset[0]*ug[4];
-
 
582
-
 
583
                ul[0] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
-
 
584
                ul[1] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
-
 
585
                ul[2] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
-
 
-
 
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];
586
        }
725
        }
587
-
 
-
 
726
       
588
        if (nodeJOffset) {
727
        if (nodeJOffset) {
589
                Wu[0] =  nodeJOffset[2]*ug[10] - nodeJOffset[1]*ug[11];
-
 
590
                Wu[1] = -nodeJOffset[2]*ug[9]  + nodeJOffset[0]*ug[11];
-
 
591
                Wu[2] =  nodeJOffset[1]*ug[9]  - nodeJOffset[0]*ug[10];
-
 
592
-
 
593
                ul[6] += R[0][0]*Wu[0] + R[0][1]*Wu[1] + R[0][2]*Wu[2];
-
 
594
                ul[7] += R[1][0]*Wu[0] + R[1][1]*Wu[1] + R[1][2]*Wu[2];
-
 
595
                ul[8] += R[2][0]*Wu[0] + R[2][1]*Wu[1] + R[2][2]*Wu[2];
-
 
-
 
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];
596
        }
735
        }
597
-
 
598
        ub(0) = ul[6] - ul[0];
-
 
-
 
736
       
-
 
737
        ab(0) = al[6] - al[0];
599
        double tmp;
738
        double tmp;
600
        tmp = oneOverL*(ul[1]-ul[7]);
-
 
601
        ub(1) = ul[5] + tmp;
-
 
602
        ub(2) = ul[11] + tmp;
-
 
603
        tmp = oneOverL*(ul[8]-ul[2]);
-
 
604
        ub(3) = ul[4] + tmp;
-
 
605
        ub(4) = ul[10] + tmp;
-
 
606
        ub(5) = ul[9] - ul[3];
-
 
607
-
 
608
        return ub;
-
 
-
 
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;
609
}
748
}
610
749
611
750
612
const Vector &
751
const Vector &
613
PDeltaCrdTransf3d::getGlobalResistingForce(const Vector &pb, const Vector &p0)
752
PDeltaCrdTransf3d::getGlobalResistingForce(const Vector &pb, const Vector &p0)
614
{
753
{
615
        // transform resisting forces from the basic system to local coordinates
-
 
616
        static double pl[12];
-
 
617
-
 
618
        double q0 = pb(0);
-
 
619
        double q1 = pb(1);
-
 
620
        double q2 = pb(2);
-
 
621
        double q3 = pb(3);
-
 
622
        double q4 = pb(4);
-
 
623
        double q5 = pb(5);
-
 
624
-
 
625
        double oneOverL = 1.0/L;
-
 
626
-
 
627
        pl[0]  = -q0;
-
 
628
        pl[1]  =  oneOverL*(q1+q2);
-
 
629
        pl[2]  = -oneOverL*(q3+q4);
-
 
630
        pl[3]  = -q5;
-
 
631
        pl[4]  =  q3;
-
 
632
        pl[5]  =  q1;
-
 
633
        pl[6]  =  q0;
-
 
634
        pl[7]  = -pl[1];
-
 
635
        pl[8]  = -pl[2];
-
 
636
        pl[9]  =  q5;
-
 
637
        pl[10] =  q4;
-
 
638
        pl[11] =  q2;
-
 
639
-
 
640
        pl[0] += p0(0);
-
 
641
        pl[1] += p0(1);
-
 
642
        pl[7] += p0(2);
-
 
643
        pl[2] += p0(3);
-
 
644
        pl[8] += p0(4);
-
 
645
-
 
646
        // Include leaning column effects (P-Delta)
-
 
647
        double NoverL;
-
 
648
        NoverL = ul17*q0*oneOverL;            
-
 
649
        pl[1] += NoverL;
-
 
650
        pl[7] -= NoverL;
-
 
651
        NoverL = ul28*q0*oneOverL;
-
 
652
        pl[2] += NoverL;
-
 
653
        pl[8] -= NoverL;
-
 
654
-
 
655
        // transform resisting forces  from local to global coordinates
-
 
656
        static Vector pg(12);
-
 
657
-
 
658
        pg(0)  = R[0][0]*pl[0] + R[1][0]*pl[1] + R[2][0]*pl[2];
-
 
659
        pg(1)  = R[0][1]*pl[0] + R[1][1]*pl[1] + R[2][1]*pl[2];
-
 
660
        pg(2)  = R[0][2]*pl[0] + R[1][2]*pl[1] + R[2][2]*pl[2];
-
 
661
-
 
662
        pg(3)  = R[0][0]*pl[3] + R[1][0]*pl[4] + R[2][0]*pl[5];
-
 
663
        pg(4)  = R[0][1]*pl[3] + R[1][1]*pl[4] + R[2][1]*pl[5];
-
 
664
        pg(5)  = R[0][2]*pl[3] + R[1][2]*pl[4] + R[2][2]*pl[5];
-
 
665
-
 
666
        pg(6)  = R[0][0]*pl[6] + R[1][0]*pl[7] + R[2][0]*pl[8];
-
 
667
        pg(7)  = R[0][1]*pl[6] + R[1][1]*pl[7] + R[2][1]*pl[8];
-
 
668
        pg(8)  = R[0][2]*pl[6] + R[1][2]*pl[7] + R[2][2]*pl[8];
-
 
669
-
 
670
        pg(9)  = R[0][0]*pl[9] + R[1][0]*pl[10] + R[2][0]*pl[11];
-
 
671
        pg(10) = R[0][1]*pl[9] + R[1][1]*pl[10] + R[2][1]*pl[11];
-
 
672
        pg(11) = R[0][2]*pl[9] + R[1][2]*pl[10] + R[2][2]*pl[11];
-
 
673
-
 
674
        if (nodeIOffset) {
-
 
675
                pg(3) += -nodeIOffset[2]*pg(1) + nodeIOffset[1]*pg(2);
-
 
676
                pg(4) +=  nodeIOffset[2]*pg(0) - nodeIOffset[0]*pg(2);
-
 
677
                pg(5) += -nodeIOffset[1]*pg(0) + nodeIOffset[0]*pg(1);
-
 
678
        }
-
 
679
-
 
680
        if (nodeJOffset) {
-
 
681
                pg(9)  += -nodeJOffset[2]*pg(7) + nodeJOffset[1]*pg(8);
-
 
682
                pg(10) +=  nodeJOffset[2]*pg(6) - nodeJOffset[0]*pg(8);
-
 
683
                pg(11) += -nodeJOffset[1]*pg(6) + nodeJOffset[0]*pg(7);
-
 
684
        }
-
 
685
       
-
 
686
        return pg;
-
 
-
 
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;
687
}
826
}
688
827
689
828
690
const Matrix &
829
const Matrix &
691
PDeltaCrdTransf3d::getGlobalStiffMatrix (const Matrix &KB, const Vector &pb)
830
PDeltaCrdTransf3d::getGlobalStiffMatrix (const Matrix &KB, const Vector &pb)
692
{
831
{
693
        static Matrix kg(12,12);        // Global stiffness for return
-
 
694
        static double kb[6][6];         // Basic stiffness
-
 
695
        static double kl[12][12];       // Local stiffness
-
 
696
        static double tmp[12][12];      // Temporary storage
-
 
697
-
 
698
        double oneOverL = 1.0/L;
-
 
699
-
 
700
        int i,j;
-
 
701
        for (i = 0; i < 6; i++)
-
 
702
                for (j = 0; j < 6; j++)
-
 
703
                        kb[i][j] = KB(i,j);
-
 
704
-
 
705
        // Transform basic stiffness to local system
-
 
706
        // First compute kb*T_{bl}
-
 
707
        for (i = 0; i < 6; i++) {
-
 
708
                tmp[i][0]  = -kb[i][0];
-
 
709
                tmp[i][1]  =  oneOverL*(kb[i][1]+kb[i][2]);
-
 
710
                tmp[i][2]  = -oneOverL*(kb[i][3]+kb[i][4]);
-
 
711
                tmp[i][3]  = -kb[i][5];
-
 
712
                tmp[i][4]  =  kb[i][3];
-
 
713
                tmp[i][5]  =  kb[i][1];
-
 
714
                tmp[i][6]  =  kb[i][0];
-
 
715
                tmp[i][7]  = -tmp[i][1];
-
 
716
                tmp[i][8]  = -tmp[i][2];
-
 
717
                tmp[i][9]  =  kb[i][5];
-
 
718
                tmp[i][10] =  kb[i][4];
-
 
719
                tmp[i][11] =  kb[i][2];
-
 
720
        }
-
 
721
       
-
 
722
        // Now compute T'_{bl}*(kb*T_{bl})
-
 
723
        for (i = 0; i < 12; i++) {
-
 
724
                kl[0][i]  = -tmp[0][i];
-
 
725
                kl[1][i]  =  oneOverL*(tmp[1][i]+tmp[2][i]);
-
 
726
                kl[2][i]  = -oneOverL*(tmp[3][i]+tmp[4][i]);
-
 
727
                kl[3][i]  = -tmp[5][i];
-
 
728
                kl[4][i]  =  tmp[3][i];
-
 
729
                kl[5][i]  =  tmp[1][i];
-
 
730
                kl[6][i]  =  tmp[0][i];
-
 
731
                kl[7][i]  = -kl[1][i];
-
 
732
                kl[8][i]  = -kl[2][i];
-
 
733
                kl[9][i]  =  tmp[5][i];
-
 
734
                kl[10][i] =  tmp[4][i];
-
 
735
                kl[11][i] =  tmp[2][i];
-
 
736
        }
-
 
737
-
 
738
        // Include geometric stiffness effects in local system
-
 
739
        double NoverL = pb(0)*oneOverL;
-
 
740
        kl[1][1] += NoverL;
-
 
741
        kl[2][2] += NoverL;
-
 
742
        kl[7][7] += NoverL;
-
 
743
        kl[8][8] += NoverL;
-
 
744
        kl[1][7] -= NoverL;
-
 
745
        kl[7][1] -= NoverL;
-
 
746
        kl[2][8] -= NoverL;
-
 
747
        kl[8][2] -= NoverL;
-
 
748
       
-
 
749
        static double RWI[3][3];
-
 
750
-
 
751
        if (nodeIOffset) {
-
 
752
                // Compute RWI
-
 
753
                RWI[0][0] = -R[0][1]*nodeIOffset[2] + R[0][2]*nodeIOffset[1];
-
 
754
                RWI[1][0] = -R[1][1]*nodeIOffset[2] + R[1][2]*nodeIOffset[1];
-
 
755
                RWI[2][0] = -R[2][1]*nodeIOffset[2] + R[2][2]*nodeIOffset[1];
-
 
756
-
 
757
                RWI[0][1] =  R[0][0]*nodeIOffset[2] - R[0][2]*nodeIOffset[0];
-
 
758
                RWI[1][1] =  R[1][0]*nodeIOffset[2] - R[1][2]*nodeIOffset[0];
-
 
759
                RWI[2][1] =  R[2][0]*nodeIOffset[2] - R[2][2]*nodeIOffset[0];
-
 
760
-
 
761
                RWI[0][2] = -R[0][0]*nodeIOffset[1] + R[0][1]*nodeIOffset[0];
-
 
762
                RWI[1][2] = -R[1][0]*nodeIOffset[1] + R[1][1]*nodeIOffset[0];
-
 
763
                RWI[2][2] = -R[2][0]*nodeIOffset[1] + R[2][1]*nodeIOffset[0];
-
 
764
        }
-
 
765
-
 
766
        static double RWJ[3][3];
-
 
767
-
 
768
        if (nodeJOffset) {
-
 
769
                // Compute RWJ
-
 
770
                RWJ[0][0] = -R[0][1]*nodeJOffset[2] + R[0][2]*nodeJOffset[1];
-
 
771
                RWJ[1][0] = -R[1][1]*nodeJOffset[2] + R[1][2]*nodeJOffset[1];
-
 
772
                RWJ[2][0] = -R[2][1]*nodeJOffset[2] + R[2][2]*nodeJOffset[1];
-
 
773
-
 
774
                RWJ[0][1] =  R[0][0]*nodeJOffset[2] - R[0][2]*nodeJOffset[0];
-
 
775
                RWJ[1][1] =  R[1][0]*nodeJOffset[2] - R[1][2]*nodeJOffset[0];
-
 
776
                RWJ[2][1] =  R[2][0]*nodeJOffset[2] - R[2][2]*nodeJOffset[0];
-
 
777
-
 
778
                RWJ[0][2] = -R[0][0]*nodeJOffset[1] + R[0][1]*nodeJOffset[0];
-
 
779
                RWJ[1][2] = -R[1][0]*nodeJOffset[1] + R[1][1]*nodeJOffset[0];
-
 
780
                RWJ[2][2] = -R[2][0]*nodeJOffset[1] + R[2][1]*nodeJOffset[0];
-
 
781
        }
-
 
782
-
 
783
        // Transform local stiffness to global system
-
 
784
        // First compute kl*T_{lg}
-
 
785
        int m;
-
 
786
        for (m = 0; m < 12; m++) {
-
 
787
                tmp[m][0] = kl[m][0]*R[0][0] + kl[m][1]*R[1][0]  + kl[m][2]*R[2][0];
-
 
788
                tmp[m][1] = kl[m][0]*R[0][1] + kl[m][1]*R[1][1]  + kl[m][2]*R[2][1];
-
 
789
                tmp[m][2] = kl[m][0]*R[0][2] + kl[m][1]*R[1][2]  + kl[m][2]*R[2][2];
-
 
790
-
 
791
                tmp[m][3] = kl[m][3]*R[0][0] + kl[m][4]*R[1][0]  + kl[m][5]*R[2][0];
-
 
792
                tmp[m][4] = kl[m][3]*R[0][1] + kl[m][4]*R[1][1]  + kl[m][5]*R[2][1];
-
 
793
                tmp[m][5] = kl[m][3]*R[0][2] + kl[m][4]*R[1][2]  + kl[m][5]*R[2][2];
-
 
794
-
 
795
                if (nodeIOffset) {
-
 
796
                        tmp[m][3]  += kl[m][0]*RWI[0][0]  + kl[m][1]*RWI[1][0]  + kl[m][2]*RWI[2][0];
-
 
797
                        tmp[m][4]  += kl[m][0]*RWI[0][1]  + kl[m][1]*RWI[1][1]  + kl[m][2]*RWI[2][1];
-
 
798
                        tmp[m][5]  += kl[m][0]*RWI[0][2]  + kl[m][1]*RWI[1][2]  + kl[m][2]*RWI[2][2];
-
 
799
                }
-
 
800
-
 
801
                tmp[m][6] = kl[m][6]*R[0][0] + kl[m][7]*R[1][0]  + kl[m][8]*R[2][0];
-
 
802
                tmp[m][7] = kl[m][6]*R[0][1] + kl[m][7]*R[1][1]  + kl[m][8]*R[2][1];
-
 
803
                tmp[m][8] = kl[m][6]*R[0][2] + kl[m][7]*R[1][2]  + kl[m][8]*R[2][2];
-
 
804
-
 
805
                tmp[m][9]  = kl[m][9]*R[0][0] + kl[m][10]*R[1][0] + kl[m][11]*R[2][0];
-
 
806
                tmp[m][10] = kl[m][9]*R[0][1] + kl[m][10]*R[1][1] + kl[m][11]*R[2][1];
-
 
807
                tmp[m][11] = kl[m][9]*R[0][2] + kl[m][10]*R[1][2] + kl[m][11]*R[2][2];
-
 
808
-
 
809
                if (nodeJOffset) {
-
 
810
                        tmp[m][9]   += kl[m][6]*RWJ[0][0]  + kl[m][7]*RWJ[1][0]  + kl[m][8]*RWJ[2][0];
-
 
811
                        tmp[m][10]  += kl[m][6]*RWJ[0][1]  + kl[m][7]*RWJ[1][1]  + kl[m][8]*RWJ[2][1];
-
 
812
                        tmp[m][11]  += kl[m][6]*RWJ[0][2]  + kl[m][7]*RWJ[1][2]  + kl[m][8]*RWJ[2][2];
-
 
813
                }
-
 
814
-
 
815
        }
-
 
816
-
 
817
        // Now compute T'_{lg}*(kl*T_{lg})
-
 
818
        for (m = 0; m < 12; m++) {
-
 
819
                kg(0,m) = R[0][0]*tmp[0][m] + R[1][0]*tmp[1][m]  + R[2][0]*tmp[2][m];
-
 
820
                kg(1,m) = R[0][1]*tmp[0][m] + R[1][1]*tmp[1][m]  + R[2][1]*tmp[2][m];
-
 
821
                kg(2,m) = R[0][2]*tmp[0][m] + R[1][2]*tmp[1][m]  + R[2][2]*tmp[2][m];
-
 
822
-
 
823
                kg(3,m) = R[0][0]*tmp[3][m] + R[1][0]*tmp[4][m]  + R[2][0]*tmp[5][m];
-
 
824
                kg(4,m) = R[0][1]*tmp[3][m] + R[1][1]*tmp[4][m]  + R[2][1]*tmp[5][m];
-
 
825
                kg(5,m) = R[0][2]*tmp[3][m] + R[1][2]*tmp[4][m]  + R[2][2]*tmp[5][m];
-
 
826
-
 
827
                if (nodeIOffset) {
-
 
828
                        kg(3,m) += RWI[0][0]*tmp[0][m]  + RWI[1][0]*tmp[1][m] + RWI[2][0]*tmp[2][m];
-
 
829
                        kg(4,m) += RWI[0][1]*tmp[0][m]  + RWI[1][1]*tmp[1][m] + RWI[2][1]*tmp[2][m];
-
 
830
                        kg(5,m) += RWI[0][2]*tmp[0][m]  + RWI[1][2]*tmp[1][m] + RWI[2][2]*tmp[2][m];
-
 
831
                }
-
 
832
-
 
833
                kg(6,m) = R[0][0]*tmp[6][m] + R[1][0]*tmp[7][m]  + R[2][0]*tmp[8][m];
-
 
834
                kg(7,m) = R[0][1]*tmp[6][m] + R[1][1]*tmp[7][m]  + R[2][1]*tmp[8][m];
-
 
835
                kg(8,m) = R[0][2]*tmp[6][m] + R[1][2]*tmp[7][m]  + R[2][2]*tmp[8][m];
-
 
836
-
 
837
                kg(9,m)  = R[0][0]*tmp[9][m] + R[1][0]*tmp[10][m] + R[2][0]*tmp[11][m];
-
 
838
                kg(10,m) = R[0][1]*tmp[9][m] + R[1][1]*tmp[10][m] + R[2][1]*tmp[11][m];
-
 
839
                kg(11,m) = R[0][2]*tmp[9][m] + R[1][2]*tmp[10][m] + R[2][2]*tmp[11][m];
-
 
840
-
 
841
                if (nodeJOffset) {
-
 
842
                        kg(9,m)  += RWJ[0][0]*tmp[6][m]  + RWJ[1][0]*tmp[7][m] + RWJ[2][0]*tmp[8][m];
-
 
843
                        kg(10,m) += RWJ[0][1]*tmp[6][m]  + RWJ[1][1]*tmp[7][m] + RWJ[2][1]*tmp[8][m];
-
 
844
                        kg(11,m) += RWJ[0][2]*tmp[6][m]  + RWJ[1][2]*tmp[7][m] + RWJ[2][2]*tmp[8][m];
-
 
845
                }
-
 
846
        }
-
 
847
-
 
848
        return kg;
-
 
-
 
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++)