WebSVN - OpenSees - Blame - Rev 755 - /trunk/SRC/coordTransformation/LinearCrdTransf3d.cpp

Rev	Author	Line No.	Line
2	fmk	1	/* ****************************************************************
		2	OpenSees - Open System for Earthquake Engineering Simulation
		3	Pacific Earthquake Engineering Research Center
		4
		5
		6	(C) Copyright 1999, The Regents of the University of California
		7	All Rights Reserved.
		8
		9	Commercial use of this program without express permission of the
		10	University of California, Berkeley, is strictly prohibited. See
		11	file 'COPYRIGHT' in main directory for information on usage and
		12	redistribution, and for a DISCLAIMER OF ALL WARRANTIES.
		13
		14	Developed by:
		15	Frank McKenna (fmckenna@ce.berkeley.edu)
		16	Gregory L. Fenves (fenves@ce.berkeley.edu)
		17	Filip C. Filippou (filippou@ce.berkeley.edu)
		18
		19	**************************************************************** */
		20
755	jeremic	21	// $Revision: 1.5 $
		22	// $Date: 2002-01-06 19:22:27 $
2	fmk	23	// $Source: /usr/local/cvs/OpenSees/SRC/coordTransformation/LinearCrdTransf3d.cpp,v $
		24
		25
272	mhscott	26	// File: ~/crdTransf/LinearCrdTransf3d.C
2	fmk	27	//
		28	// Written: Remo Magalhaes de Souza (rmsouza@ce.berkeley.edu)
		29	// Created: 04/2000
		30	// Revision: A
		31	//
		32	// Purpose: This file contains the implementation for the
		33	// LinearCrdTransf3d class. LinearCrdTransf3d is a linear
272	mhscott	34	// transformation for a planar frame between the global
2	fmk	35	// and basic coordinate systems
		36
		37
		38	#include <Vector.h>
		39	#include <Matrix.h>
		40	#include <Node.h>
		41	#include <Channel.h>
		42
272	mhscott	43	#include <iomanip.h>
		44
2	fmk	45	#include <LinearCrdTransf3d.h>
		46
272	mhscott	47	// constructor:
		48	LinearCrdTransf3d::LinearCrdTransf3d(int tag, const Vector &vecInLocXZPlane):
		49	CrdTransf3d(tag, CRDTR_TAG_LinearCrdTransf3d),
		50	nodeIPtr(0), nodeJPtr(0),
		51	nodeIOffset(0), nodeJOffset(0), L(0)
		52	{
610	mhscott	53	for (int i = 0; i < 2; i++)
		54	for (int j = 0; j < 3; j++)
		55	R[i][j] = 0.0;
2	fmk	56
610	mhscott	57	R[2][0] = vecInLocXZPlane(0);
		58	R[2][1] = vecInLocXZPlane(1);
		59	R[2][2] = vecInLocXZPlane(2);
		60
272	mhscott	61	// Does nothing
		62	}
		63
2	fmk	64	// constructor:
		65	LinearCrdTransf3d::LinearCrdTransf3d(int tag, const Vector &vecInLocXZPlane,
272	mhscott	66	const Vector &rigJntOffset1,
		67	const Vector &rigJntOffset2):
2	fmk	68	CrdTransf3d(tag, CRDTR_TAG_LinearCrdTransf3d),
		69	nodeIPtr(0), nodeJPtr(0),
272	mhscott	70	nodeIOffset(0), nodeJOffset(0), L(0)
2	fmk	71	{
610	mhscott	72	for (int i = 0; i < 2; i++)
		73	for (int j = 0; j < 3; j++)
		74	R[i][j] = 0.0;
272	mhscott	75
610	mhscott	76	R[2][0] = vecInLocXZPlane(0);
		77	R[2][1] = vecInLocXZPlane(1);
		78	R[2][2] = vecInLocXZPlane(2);
		79
272	mhscott	80	// check rigid joint offset for node I
		81	if (&rigJntOffset1 == 0 \|\| rigJntOffset1.Size() != 3 ) {
		82	cerr << "LinearCrdTransf3d::LinearCrdTransf3d: Invalid rigid joint offset vector for node I\n";
		83	cerr << "Size must be 3\n";
		84	}
		85	else if (rigJntOffset1.Norm() > 0.0) {
		86	nodeIOffset = new double[3];
		87	nodeIOffset[0] = rigJntOffset1(0);
		88	nodeIOffset[1] = rigJntOffset1(1);
		89	nodeIOffset[2] = rigJntOffset1(2);
		90	}
2	fmk	91
		92	// check rigid joint offset for node J
272	mhscott	93	if (&rigJntOffset2 == 0 \|\| rigJntOffset2.Size() != 3 ) {
		94	cerr << "LinearCrdTransf3d::LinearCrdTransf3d: Invalid rigid joint offset vector for node J\n";
		95	cerr << "Size must be 3\n";
		96	}
		97	else if (rigJntOffset2.Norm() > 0.0) {
		98	nodeJOffset = new double[3];
		99	nodeJOffset[0] = rigJntOffset2(0);
		100	nodeJOffset[1] = rigJntOffset2(1);
		101	nodeJOffset[2] = rigJntOffset2(2);
		102	}
2	fmk	103	}
		104
		105
		106
272	mhscott	107
2	fmk	108	// constructor:
		109	// invoked by a FEM_ObjectBroker, recvSelf() needs to be invoked on this object.
		110	LinearCrdTransf3d::LinearCrdTransf3d():
		111	CrdTransf3d(0, CRDTR_TAG_LinearCrdTransf3d),
		112	nodeIPtr(0), nodeJPtr(0),
272	mhscott	113	nodeIOffset(0), nodeJOffset(0), L(0)
2	fmk	114	{
610	mhscott	115	for (int i = 0; i < 3; i++)
		116	for (int j = 0; j < 3; j++)
		117	R[i][j] = 0.0;
2	fmk	118	}
		119
		120
		121
		122	// destructor:
		123	LinearCrdTransf3d::~LinearCrdTransf3d()
		124	{
272	mhscott	125	if (nodeIOffset)
		126	delete [] nodeIOffset;
		127	if (nodeJOffset)
		128	delete [] nodeJOffset;
2	fmk	129	}
		130
272	mhscott	131
2	fmk	132	int
		133	LinearCrdTransf3d::commitState(void)
		134	{
		135	return 0;
		136	}
		137
		138
		139	int
		140	LinearCrdTransf3d::revertToLastCommit(void)
		141	{
		142	return 0;
		143	}
		144
		145
		146	int
		147	LinearCrdTransf3d::revertToStart(void)
		148	{
		149	return 0;
		150	}
		151
		152
272	mhscott	153	int
		154	LinearCrdTransf3d::initialize(Node nodeIPointer, Node nodeJPointer)
2	fmk	155	{
		156	int error;
		157
		158	nodeIPtr = nodeIPointer;
		159	nodeJPtr = nodeJPointer;
		160
		161	if ((!nodeIPtr) \|\| (!nodeJPtr))
		162	{
		163	cerr << "\nLinearCrdTransf3d::initialize";
		164	cerr << "\ninvalid pointers to the element nodes\n";
		165	return -1;
		166	}
		167
		168	// get element length and orientation
		169	if ((error = this->computeElemtLengthAndOrient()))
		170	return error;
		171
272	mhscott	172	// get 3by3 rotation matrix
2	fmk	173	if ((error = this->getLocalAxes()))
		174	return error;
		175
		176	return 0;
		177	}
		178
		179
		180	int
		181	LinearCrdTransf3d::update(void)
272	mhscott	182	{
		183	return 0;
2	fmk	184	}
		185
		186
		187	int
		188	LinearCrdTransf3d::computeElemtLengthAndOrient()
		189	{
		190	// element projection
		191	static Vector dx(3);
		192
272	mhscott	193	const Vector &ndICoords = nodeIPtr->getCrds();
		194	const Vector &ndJCoords = nodeJPtr->getCrds();
		195
		196	if (nodeIOffset == 0) {
		197	dx(0) = ndJCoords(0) - ndICoords(0);
		198	dx(1) = ndJCoords(1) - ndICoords(1);
		199	dx(2) = ndJCoords(2) - ndICoords(2);
		200	}
		201	else {
		202	dx(0) = ndJCoords(0) + nodeJOffset[0] - ndICoords(0) - nodeIOffset[0];
		203	dx(1) = ndJCoords(1) + nodeJOffset[1] - ndICoords(1) - nodeIOffset[1];
		204	dx(2) = ndJCoords(2) + nodeJOffset[2] - ndICoords(2) - nodeIOffset[2];
		205	}
		206
2	fmk	207	// calculate the element length
		208	L = dx.Norm();
		209
272	mhscott	210	if (L == 0.0) {
2	fmk	211	cerr << "\nLinearCrdTransf3d::computeElemtLengthAndOrien: 0 length\n";
		212	return -2;
		213	}
		214
		215	// calculate the element local x axis components (direction cossines)
272	mhscott	216	// wrt to the global coordinates
		217	R[0][0] = dx(0)/L;
		218	R[0][1] = dx(1)/L;
		219	R[0][2] = dx(2)/L;
		220
2	fmk	221	return 0;
		222	}
		223
		224
272	mhscott	225	int
		226	LinearCrdTransf3d::getLocalAxes(void)
2	fmk	227	{
272	mhscott	228	// Compute y = v cross x
		229	// Note: v(i) is stored in R[2][i]
		230	static Vector vAxis(3);
		231	vAxis(0) = R[2][0]; vAxis(1) = R[2][1]; vAxis(2) = R[2][2];
2	fmk	232
272	mhscott	233	static Vector xAxis(3);
		234	xAxis(0) = R[0][0]; xAxis(1) = R[0][1]; xAxis(2) = R[0][2];
2	fmk	235
272	mhscott	236	static Vector yAxis(3);
2	fmk	237
272	mhscott	238	yAxis(0) = vAxis(1)xAxis(2) - vAxis(2)xAxis(1);
		239	yAxis(1) = vAxis(2)xAxis(0) - vAxis(0)xAxis(2);
		240	yAxis(2) = vAxis(0)xAxis(1) - vAxis(1)xAxis(0);
2	fmk	241
272	mhscott	242	double ynorm = yAxis.Norm();
		243
		244	if (ynorm == 0) {
		245	cerr << "\nLinearCrdTransf3d::getLocalAxes";
		246	cerr << "\nvector v that defines plane xz is parallel to x axis\n";
		247	return -3;
		248	}
		249
		250	yAxis /= ynorm;
		251
		252	// Compute z = x cross y
		253	static Vector zAxis(3);
		254
		255	zAxis(0) = xAxis(1)yAxis(2) - xAxis(2)yAxis(1);
		256	zAxis(1) = xAxis(2)yAxis(0) - xAxis(0)yAxis(2);
		257	zAxis(2) = xAxis(0)yAxis(1) - xAxis(1)yAxis(0);
		258
		259	// Fill in transformation matrix
		260	R[1][0] = yAxis(0);
		261	R[1][1] = yAxis(1);
		262	R[1][2] = yAxis(2);
		263
		264	R[2][0] = zAxis(0);
		265	R[2][1] = zAxis(1);
		266	R[2][2] = zAxis(2);
		267
		268	return 0;
2	fmk	269	}
		270
		271
272	mhscott	272
2	fmk	273	double
		274	LinearCrdTransf3d::getInitialLength(void)
		275	{
		276	return L;
		277	}
		278
		279
		280	double
		281	LinearCrdTransf3d::getDeformedLength(void)
		282	{
		283	return L;
		284	}
		285
		286
		287	const Vector &
		288	LinearCrdTransf3d::getBasicTrialDisp (void)
		289	{
272	mhscott	290	// determine global displacements
		291	const Vector &disp1 = nodeIPtr->getTrialDisp();
		292	const Vector &disp2 = nodeJPtr->getTrialDisp();
2	fmk	293
272	mhscott	294	static double ug[12];
		295	for (int i = 0; i < 6; i++) {
		296	ug[i] = disp1(i);
		297	ug[i+6] = disp2(i);
		298	}
2	fmk	299
272	mhscott	300	double oneOverL = 1.0/L;
2	fmk	301
272	mhscott	302	static Vector ub(6);
2	fmk	303
272	mhscott	304	static double ul[12];
2	fmk	305
272	mhscott	306	ul[0] = R[0][0]ug[0] + R[0][1]ug[1] + R[0][2]*ug[2];
		307	ul[1] = R[1][0]ug[0] + R[1][1]ug[1] + R[1][2]*ug[2];
		308	ul[2] = R[2][0]ug[0] + R[2][1]ug[1] + R[2][2]*ug[2];
		309
		310	ul[3] = R[0][0]ug[3] + R[0][1]ug[4] + R[0][2]*ug[5];
		311	ul[4] = R[1][0]ug[3] + R[1][1]ug[4] + R[1][2]*ug[5];
		312	ul[5] = R[2][0]ug[3] + R[2][1]ug[4] + R[2][2]*ug[5];
		313
		314	ul[6] = R[0][0]ug[6] + R[0][1]ug[7] + R[0][2]*ug[8];
		315	ul[7] = R[1][0]ug[6] + R[1][1]ug[7] + R[1][2]*ug[8];
		316	ul[8] = R[2][0]ug[6] + R[2][1]ug[7] + R[2][2]*ug[8];
		317
		318	ul[9] = R[0][0]ug[9] + R[0][1]ug[10] + R[0][2]*ug[11];
		319	ul[10] = R[1][0]ug[9] + R[1][1]ug[10] + R[1][2]*ug[11];
		320	ul[11] = R[2][0]ug[9] + R[2][1]ug[10] + R[2][2]*ug[11];
		321
		322	static double Wu[3];
		323	if (nodeIOffset) {
		324	Wu[0] = nodeIOffset[2]ug[4] - nodeIOffset[1]ug[5];
		325	Wu[1] = -nodeIOffset[2]ug[3] + nodeIOffset[0]ug[5];
		326	Wu[2] = nodeIOffset[1]ug[3] - nodeIOffset[0]ug[4];
		327
		328	ul[0] += R[0][0]Wu[0] + R[0][1]Wu[1] + R[0][2]*Wu[2];
		329	ul[1] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		330	ul[2] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		331	}
		332
		333	if (nodeJOffset) {
		334	Wu[0] = nodeJOffset[2]ug[10] - nodeJOffset[1]ug[11];
		335	Wu[1] = -nodeJOffset[2]ug[9] + nodeJOffset[0]ug[11];
		336	Wu[2] = nodeJOffset[1]ug[9] - nodeJOffset[0]ug[10];
		337
		338	ul[6] += R[0][0]Wu[0] + R[0][1]Wu[1] + R[0][2]*Wu[2];
		339	ul[7] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		340	ul[8] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		341	}
		342
		343	ub(0) = ul[6] - ul[0];
		344	double tmp;
		345	tmp = oneOverL*(ul[1]-ul[7]);
		346	ub(1) = ul[5] + tmp;
		347	ub(2) = ul[11] + tmp;
		348	tmp = oneOverL*(ul[8]-ul[2]);
		349	ub(3) = ul[4] + tmp;
		350	ub(4) = ul[10] + tmp;
		351	ub(5) = ul[9] - ul[3];
		352
		353	return ub;
2	fmk	354	}
		355
		356
		357	const Vector &
		358	LinearCrdTransf3d::getBasicIncrDisp (void)
		359	{
272	mhscott	360	// determine global displacements
		361	const Vector &disp1 = nodeIPtr->getIncrDisp();
		362	const Vector &disp2 = nodeJPtr->getIncrDisp();
2	fmk	363
272	mhscott	364	static double ug[12];
		365	for (int i = 0; i < 6; i++) {
		366	ug[i] = disp1(i);
		367	ug[i+6] = disp2(i);
		368	}
2	fmk	369
272	mhscott	370	double oneOverL = 1.0/L;
2	fmk	371
272	mhscott	372	static Vector ub(6);
2	fmk	373
272	mhscott	374	static double ul[12];
2	fmk	375
272	mhscott	376	ul[0] = R[0][0]ug[0] + R[0][1]ug[1] + R[0][2]*ug[2];
		377	ul[1] = R[1][0]ug[0] + R[1][1]ug[1] + R[1][2]*ug[2];
		378	ul[2] = R[2][0]ug[0] + R[2][1]ug[1] + R[2][2]*ug[2];
		379
		380	ul[3] = R[0][0]ug[3] + R[0][1]ug[4] + R[0][2]*ug[5];
		381	ul[4] = R[1][0]ug[3] + R[1][1]ug[4] + R[1][2]*ug[5];
		382	ul[5] = R[2][0]ug[3] + R[2][1]ug[4] + R[2][2]*ug[5];
		383
		384	ul[6] = R[0][0]ug[6] + R[0][1]ug[7] + R[0][2]*ug[8];
		385	ul[7] = R[1][0]ug[6] + R[1][1]ug[7] + R[1][2]*ug[8];
		386	ul[8] = R[2][0]ug[6] + R[2][1]ug[7] + R[2][2]*ug[8];
		387
		388	ul[9] = R[0][0]ug[9] + R[0][1]ug[10] + R[0][2]*ug[11];
		389	ul[10] = R[1][0]ug[9] + R[1][1]ug[10] + R[1][2]*ug[11];
		390	ul[11] = R[2][0]ug[9] + R[2][1]ug[10] + R[2][2]*ug[11];
		391
		392	static double Wu[3];
		393	if (nodeIOffset) {
		394	Wu[0] = nodeIOffset[2]ug[4] - nodeIOffset[1]ug[5];
		395	Wu[1] = -nodeIOffset[2]ug[3] + nodeIOffset[0]ug[5];
		396	Wu[2] = nodeIOffset[1]ug[3] - nodeIOffset[0]ug[4];
		397
		398	ul[0] += R[0][0]Wu[0] + R[0][1]Wu[1] + R[0][2]*Wu[2];
		399	ul[1] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		400	ul[2] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		401	}
		402
		403	if (nodeJOffset) {
		404	Wu[0] = nodeJOffset[2]ug[10] - nodeJOffset[1]ug[11];
		405	Wu[1] = -nodeJOffset[2]ug[9] + nodeJOffset[0]ug[11];
		406	Wu[2] = nodeJOffset[1]ug[9] - nodeJOffset[0]ug[10];
		407
		408	ul[6] += R[0][0]Wu[0] + R[0][1]Wu[1] + R[0][2]*Wu[2];
		409	ul[7] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		410	ul[8] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		411	}
		412
		413	ub(0) = ul[6] - ul[0];
		414	double tmp;
		415	tmp = oneOverL*(ul[1]-ul[7]);
		416	ub(1) = ul[5] + tmp;
		417	ub(2) = ul[11] + tmp;
		418	tmp = oneOverL*(ul[8]-ul[2]);
		419	ub(3) = ul[4] + tmp;
		420	ub(4) = ul[10] + tmp;
		421	ub(5) = ul[9] - ul[3];
		422
		423	return ub;
2	fmk	424	}
		425
		426
		427	const Vector &
		428	LinearCrdTransf3d::getBasicIncrDeltaDisp(void)
		429	{
272	mhscott	430	// determine global displacements
		431	const Vector &disp1 = nodeIPtr->getIncrDeltaDisp();
		432	const Vector &disp2 = nodeJPtr->getIncrDeltaDisp();
2	fmk	433
272	mhscott	434	static double ug[12];
		435	for (int i = 0; i < 6; i++) {
		436	ug[i] = disp1(i);
		437	ug[i+6] = disp2(i);
		438	}
2	fmk	439
272	mhscott	440	double oneOverL = 1.0/L;
2	fmk	441
272	mhscott	442	static Vector ub(6);
2	fmk	443
272	mhscott	444	static double ul[12];
2	fmk	445
272	mhscott	446	ul[0] = R[0][0]ug[0] + R[0][1]ug[1] + R[0][2]*ug[2];
		447	ul[1] = R[1][0]ug[0] + R[1][1]ug[1] + R[1][2]*ug[2];
		448	ul[2] = R[2][0]ug[0] + R[2][1]ug[1] + R[2][2]*ug[2];
2	fmk	449
272	mhscott	450	ul[3] = R[0][0]ug[3] + R[0][1]ug[4] + R[0][2]*ug[5];
		451	ul[4] = R[1][0]ug[3] + R[1][1]ug[4] + R[1][2]*ug[5];
		452	ul[5] = R[2][0]ug[3] + R[2][1]ug[4] + R[2][2]*ug[5];
2	fmk	453
272	mhscott	454	ul[6] = R[0][0]ug[6] + R[0][1]ug[7] + R[0][2]*ug[8];
		455	ul[7] = R[1][0]ug[6] + R[1][1]ug[7] + R[1][2]*ug[8];
		456	ul[8] = R[2][0]ug[6] + R[2][1]ug[7] + R[2][2]*ug[8];
2	fmk	457
272	mhscott	458	ul[9] = R[0][0]ug[9] + R[0][1]ug[10] + R[0][2]*ug[11];
		459	ul[10] = R[1][0]ug[9] + R[1][1]ug[10] + R[1][2]*ug[11];
		460	ul[11] = R[2][0]ug[9] + R[2][1]ug[10] + R[2][2]*ug[11];
2	fmk	461
272	mhscott	462	static double Wu[3];
		463	if (nodeIOffset) {
		464	Wu[0] = nodeIOffset[2]ug[4] - nodeIOffset[1]ug[5];
		465	Wu[1] = -nodeIOffset[2]ug[3] + nodeIOffset[0]ug[5];
		466	Wu[2] = nodeIOffset[1]ug[3] - nodeIOffset[0]ug[4];
2	fmk	467
272	mhscott	468	ul[0] += R[0][0]Wu[0] + R[0][1]Wu[1] + R[0][2]*Wu[2];
		469	ul[1] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		470	ul[2] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		471	}
2	fmk	472
272	mhscott	473	if (nodeJOffset) {
		474	Wu[0] = nodeJOffset[2]ug[10] - nodeJOffset[1]ug[11];
		475	Wu[1] = -nodeJOffset[2]ug[9] + nodeJOffset[0]ug[11];
		476	Wu[2] = nodeJOffset[1]ug[9] - nodeJOffset[0]ug[10];
2	fmk	477
272	mhscott	478	ul[6] += R[0][0]Wu[0] + R[0][1]Wu[1] + R[0][2]*Wu[2];
		479	ul[7] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		480	ul[8] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		481	}
2	fmk	482
272	mhscott	483	ub(0) = ul[6] - ul[0];
		484	double tmp;
		485	tmp = oneOverL*(ul[1]-ul[7]);
		486	ub(1) = ul[5] + tmp;
		487	ub(2) = ul[11] + tmp;
		488	tmp = oneOverL*(ul[8]-ul[2]);
		489	ub(3) = ul[4] + tmp;
		490	ub(4) = ul[10] + tmp;
		491	ub(5) = ul[9] - ul[3];
2	fmk	492
272	mhscott	493	return ub;
2	fmk	494	}
		495
		496
		497	const Vector &
		498	LinearCrdTransf3d::getGlobalResistingForce(const Vector &pb, const Vector &unifLoad)
		499	{
272	mhscott	500	// transform resisting forces from the basic system to local coordinates
		501	static double pl[12];
2	fmk	502
272	mhscott	503	double q0 = pb(0);
		504	double q1 = pb(1);
		505	double q2 = pb(2);
		506	double q3 = pb(3);
		507	double q4 = pb(4);
		508	double q5 = pb(5);
2	fmk	509
272	mhscott	510	double oneOverL = 1.0/L;
		511
		512	pl[0] = -q0;
		513	pl[1] = oneOverL*(q1+q2);
		514	pl[2] = -oneOverL*(q3+q4);
		515	pl[3] = -q5;
		516	pl[4] = q3;
		517	pl[5] = q1;
		518	pl[6] = q0;
		519	pl[7] = -pl[1];
		520	pl[8] = -pl[2];
		521	pl[9] = q5;
		522	pl[10] = q4;
		523	pl[11] = q2;
		524
		525	// add end forces due to uniform distributed
		526	// loads to the system with rigid body modes
		527	pl[0] -= unifLoad(0)*L;
		528	double V;
		529	V = 0.5unifLoad(1)L;
		530	pl[1] -= V;
		531	pl[7] -= V;
		532	V = 0.5unifLoad(2)L;
		533	pl[2] -= V;
		534	pl[8] -= V;
		535
		536	// transform resisting forces from local to global coordinates
		537	static Vector pg(12);
		538
		539	pg(0) = R[0][0]pl[0] + R[1][0]pl[1] + R[2][0]*pl[2];
		540	pg(1) = R[0][1]pl[0] + R[1][1]pl[1] + R[2][1]*pl[2];
		541	pg(2) = R[0][2]pl[0] + R[1][2]pl[1] + R[2][2]*pl[2];
		542
		543	pg(3) = R[0][0]pl[3] + R[1][0]pl[4] + R[2][0]*pl[5];
		544	pg(4) = R[0][1]pl[3] + R[1][1]pl[4] + R[2][1]*pl[5];
		545	pg(5) = R[0][2]pl[3] + R[1][2]pl[4] + R[2][2]*pl[5];
		546
		547	pg(6) = R[0][0]pl[6] + R[1][0]pl[7] + R[2][0]*pl[8];
		548	pg(7) = R[0][1]pl[6] + R[1][1]pl[7] + R[2][1]*pl[8];
		549	pg(8) = R[0][2]pl[6] + R[1][2]pl[7] + R[2][2]*pl[8];
		550
		551	pg(9) = R[0][0]pl[9] + R[1][0]pl[10] + R[2][0]*pl[11];
		552	pg(10) = R[0][1]pl[9] + R[1][1]pl[10] + R[2][1]*pl[11];
		553	pg(11) = R[0][2]pl[9] + R[1][2]pl[10] + R[2][2]*pl[11];
		554
		555	if (nodeIOffset) {
		556	pg(3) += -nodeIOffset[2]pg(1) + nodeIOffset[1]pg(2);
		557	pg(4) += nodeIOffset[2]pg(0) - nodeIOffset[0]pg(2);
		558	pg(5) += -nodeIOffset[1]pg(0) + nodeIOffset[0]pg(1);
		559	}
		560
		561	if (nodeJOffset) {
		562	pg(9) += -nodeJOffset[2]pg(7) + nodeJOffset[1]pg(8);
		563	pg(10) += nodeJOffset[2]pg(6) - nodeJOffset[0]pg(8);
		564	pg(11) += -nodeJOffset[1]pg(6) + nodeJOffset[0]pg(7);
		565	}
		566
		567	return pg;
2	fmk	568	}
		569
		570	const Matrix &
272	mhscott	571	LinearCrdTransf3d::getGlobalStiffMatrix (const Matrix &KB, const Vector &pb)
2	fmk	572	{
272	mhscott	573	static Matrix kg(12,12); // Global stiffness for return
		574	static double kb[6][6]; // Basic stiffness
		575	static double kl[12][12]; // Local stiffness
		576	static double tmp[12][12]; // Temporary storage
2	fmk	577
272	mhscott	578	double oneOverL = 1.0/L;
2	fmk	579
272	mhscott	580	int i,j;
		581	for (i = 0; i < 6; i++)
		582	for (j = 0; j < 6; j++)
		583	kb[i][j] = KB(i,j);
2	fmk	584
272	mhscott	585	// Transform basic stiffness to local system
		586	// First compute kb*T_{bl}
		587	for (i = 0; i < 6; i++) {
		588	tmp[i][0] = -kb[i][0];
		589	tmp[i][1] = oneOverL*(kb[i][1]+kb[i][2]);
		590	tmp[i][2] = -oneOverL*(kb[i][3]+kb[i][4]);
		591	tmp[i][3] = -kb[i][5];
		592	tmp[i][4] = kb[i][3];
		593	tmp[i][5] = kb[i][1];
		594	tmp[i][6] = kb[i][0];
		595	tmp[i][7] = -tmp[i][1];
		596	tmp[i][8] = -tmp[i][2];
		597	tmp[i][9] = kb[i][5];
		598	tmp[i][10] = kb[i][4];
		599	tmp[i][11] = kb[i][2];
		600	}
		601
		602	// Now compute T'_{bl}(kbT_{bl})
		603	for (i = 0; i < 12; i++) {
		604	kl[0][i] = -tmp[0][i];
		605	kl[1][i] = oneOverL*(tmp[1][i]+tmp[2][i]);
		606	kl[2][i] = -oneOverL*(tmp[3][i]+tmp[4][i]);
		607	kl[3][i] = -tmp[5][i];
		608	kl[4][i] = tmp[3][i];
		609	kl[5][i] = tmp[1][i];
		610	kl[6][i] = tmp[0][i];
		611	kl[7][i] = -kl[1][i];
		612	kl[8][i] = -kl[2][i];
		613	kl[9][i] = tmp[5][i];
		614	kl[10][i] = tmp[4][i];
		615	kl[11][i] = tmp[2][i];
		616	}
2	fmk	617
272	mhscott	618	static double RWI[3][3];
		619
		620	if (nodeIOffset) {
		621	// Compute RWI
		622	RWI[0][0] = -R[0][1]nodeIOffset[2] + R[0][2]nodeIOffset[1];
		623	RWI[1][0] = -R[1][1]nodeIOffset[2] + R[1][2]nodeIOffset[1];
		624	RWI[2][0] = -R[2][1]nodeIOffset[2] + R[2][2]nodeIOffset[1];
		625
		626	RWI[0][1] = R[0][0]nodeIOffset[2] - R[0][2]nodeIOffset[0];
		627	RWI[1][1] = R[1][0]nodeIOffset[2] - R[1][2]nodeIOffset[0];
		628	RWI[2][1] = R[2][0]nodeIOffset[2] - R[2][2]nodeIOffset[0];
		629
		630	RWI[0][2] = -R[0][0]nodeIOffset[1] + R[0][1]nodeIOffset[0];
		631	RWI[1][2] = -R[1][0]nodeIOffset[1] + R[1][1]nodeIOffset[0];
		632	RWI[2][2] = -R[2][0]nodeIOffset[1] + R[2][1]nodeIOffset[0];
		633	}
		634
		635	static double RWJ[3][3];
		636
		637	if (nodeJOffset) {
		638	// Compute RWJ
		639	RWJ[0][0] = -R[0][1]nodeJOffset[2] + R[0][2]nodeJOffset[1];
		640	RWJ[1][0] = -R[1][1]nodeJOffset[2] + R[1][2]nodeJOffset[1];
		641	RWJ[2][0] = -R[2][1]nodeJOffset[2] + R[2][2]nodeJOffset[1];
		642
		643	RWJ[0][1] = R[0][0]nodeJOffset[2] - R[0][2]nodeJOffset[0];
		644	RWJ[1][1] = R[1][0]nodeJOffset[2] - R[1][2]nodeJOffset[0];
		645	RWJ[2][1] = R[2][0]nodeJOffset[2] - R[2][2]nodeJOffset[0];
		646
		647	RWJ[0][2] = -R[0][0]nodeJOffset[1] + R[0][1]nodeJOffset[0];
		648	RWJ[1][2] = -R[1][0]nodeJOffset[1] + R[1][1]nodeJOffset[0];
		649	RWJ[2][2] = -R[2][0]nodeJOffset[1] + R[2][1]nodeJOffset[0];
		650	}
		651
		652	// Transform local stiffness to global system
		653	// First compute kl*T_{lg}
		654	int m;
		655	for (m = 0; m < 12; m++) {
		656	tmp[m][0] = kl[m][0]R[0][0] + kl[m][1]R[1][0] + kl[m][2]*R[2][0];
		657	tmp[m][1] = kl[m][0]R[0][1] + kl[m][1]R[1][1] + kl[m][2]*R[2][1];
		658	tmp[m][2] = kl[m][0]R[0][2] + kl[m][1]R[1][2] + kl[m][2]*R[2][2];
		659
		660	tmp[m][3] = kl[m][3]R[0][0] + kl[m][4]R[1][0] + kl[m][5]*R[2][0];
		661	tmp[m][4] = kl[m][3]R[0][1] + kl[m][4]R[1][1] + kl[m][5]*R[2][1];
		662	tmp[m][5] = kl[m][3]R[0][2] + kl[m][4]R[1][2] + kl[m][5]*R[2][2];
		663
		664	if (nodeIOffset) {
		665	tmp[m][3] += kl[m][0]RWI[0][0] + kl[m][1]RWI[1][0] + kl[m][2]*RWI[2][0];
		666	tmp[m][4] += kl[m][0]RWI[0][1] + kl[m][1]RWI[1][1] + kl[m][2]*RWI[2][1];
		667	tmp[m][5] += kl[m][0]RWI[0][2] + kl[m][1]RWI[1][2] + kl[m][2]*RWI[2][2];
		668	}
		669
		670	tmp[m][6] = kl[m][6]R[0][0] + kl[m][7]R[1][0] + kl[m][8]*R[2][0];
		671	tmp[m][7] = kl[m][6]R[0][1] + kl[m][7]R[1][1] + kl[m][8]*R[2][1];
		672	tmp[m][8] = kl[m][6]R[0][2] + kl[m][7]R[1][2] + kl[m][8]*R[2][2];
		673
		674	tmp[m][9] = kl[m][9]R[0][0] + kl[m][10]R[1][0] + kl[m][11]*R[2][0];
		675	tmp[m][10] = kl[m][9]R[0][1] + kl[m][10]R[1][1] + kl[m][11]*R[2][1];
		676	tmp[m][11] = kl[m][9]R[0][2] + kl[m][10]R[1][2] + kl[m][11]*R[2][2];
		677
		678	if (nodeJOffset) {
		679	tmp[m][9] += kl[m][6]RWJ[0][0] + kl[m][7]RWJ[1][0] + kl[m][8]*RWJ[2][0];
		680	tmp[m][10] += kl[m][6]RWJ[0][1] + kl[m][7]RWJ[1][1] + kl[m][8]*RWJ[2][1];
		681	tmp[m][11] += kl[m][6]RWJ[0][2] + kl[m][7]RWJ[1][2] + kl[m][8]*RWJ[2][2];
		682	}
		683
		684	}
		685
		686	// Now compute T'_{lg}(klT_{lg})
		687	for (m = 0; m < 12; m++) {
		688	kg(0,m) = R[0][0]tmp[0][m] + R[1][0]tmp[1][m] + R[2][0]*tmp[2][m];
		689	kg(1,m) = R[0][1]tmp[0][m] + R[1][1]tmp[1][m] + R[2][1]*tmp[2][m];
		690	kg(2,m) = R[0][2]tmp[0][m] + R[1][2]tmp[1][m] + R[2][2]*tmp[2][m];
		691
		692	kg(3,m) = R[0][0]tmp[3][m] + R[1][0]tmp[4][m] + R[2][0]*tmp[5][m];
		693	kg(4,m) = R[0][1]tmp[3][m] + R[1][1]tmp[4][m] + R[2][1]*tmp[5][m];
		694	kg(5,m) = R[0][2]tmp[3][m] + R[1][2]tmp[4][m] + R[2][2]*tmp[5][m];
		695
		696	if (nodeIOffset) {
		697	kg(3,m) += RWI[0][0]tmp[0][m] + RWI[1][0]tmp[1][m] + RWI[2][0]*tmp[2][m];
		698	kg(4,m) += RWI[0][1]tmp[0][m] + RWI[1][1]tmp[1][m] + RWI[2][1]*tmp[2][m];
		699	kg(5,m) += RWI[0][2]tmp[0][m] + RWI[1][2]tmp[1][m] + RWI[2][2]*tmp[2][m];
		700	}
		701
		702	kg(6,m) = R[0][0]tmp[6][m] + R[1][0]tmp[7][m] + R[2][0]*tmp[8][m];
		703	kg(7,m) = R[0][1]tmp[6][m] + R[1][1]tmp[7][m] + R[2][1]*tmp[8][m];
		704	kg(8,m) = R[0][2]tmp[6][m] + R[1][2]tmp[7][m] + R[2][2]*tmp[8][m];
		705
		706	kg(9,m) = R[0][0]tmp[9][m] + R[1][0]tmp[10][m] + R[2][0]*tmp[11][m];
		707	kg(10,m) = R[0][1]tmp[9][m] + R[1][1]tmp[10][m] + R[2][1]*tmp[11][m];
		708	kg(11,m) = R[0][2]tmp[9][m] + R[1][2]tmp[10][m] + R[2][2]*tmp[11][m];
		709
		710	if (nodeJOffset) {
		711	kg(9,m) += RWJ[0][0]tmp[6][m] + RWJ[1][0]tmp[7][m] + RWJ[2][0]*tmp[8][m];
		712	kg(10,m) += RWJ[0][1]tmp[6][m] + RWJ[1][1]tmp[7][m] + RWJ[2][1]*tmp[8][m];
		713	kg(11,m) += RWJ[0][2]tmp[6][m] + RWJ[1][2]tmp[7][m] + RWJ[2][2]*tmp[8][m];
		714	}
		715	}
		716
		717	return kg;
2	fmk	718	}
272	mhscott	719
2	fmk	720
		721
		722
		723	CrdTransf3d *
		724	LinearCrdTransf3d::getCopy(void)
		725	{
		726	// create a new instance of LinearCrdTransf3d
		727
272	mhscott	728	LinearCrdTransf3d *theCopy;
		729
		730	static Vector xz(3);
		731	xz(0) = R[2][0];
		732	xz(1) = R[2][1];
		733	xz(2) = R[2][2];
		734
		735	Vector offsetI(3);
		736	Vector offsetJ(3);
		737
		738	if (nodeIOffset) {
		739	offsetI(0) = nodeIOffset[0];
		740	offsetI(1) = nodeIOffset[1];
		741	offsetI(2) = nodeIOffset[2];
		742	}
		743
		744	if (nodeJOffset) {
		745	offsetJ(0) = nodeJOffset[0];
		746	offsetJ(1) = nodeJOffset[1];
		747	offsetJ(2) = nodeJOffset[2];
		748	}
2	fmk	749
272	mhscott	750	theCopy = new LinearCrdTransf3d(this->getTag(), xz, offsetI, offsetJ);
		751
2	fmk	752	theCopy->nodeIPtr = nodeIPtr;
		753	theCopy->nodeJPtr = nodeJPtr;
		754	theCopy->L = L;
272	mhscott	755	for (int i = 0; i < 3; i++)
		756	for (int j = 0; j < 3; j++)
		757	theCopy->R[i][j] = R[i][j];
2	fmk	758
		759	return theCopy;
		760	}
		761
		762
		763	int
		764	LinearCrdTransf3d::sendSelf(int cTag, Channel &theChannel)
		765	{
		766	int res = 0;
		767
272	mhscott	768	static Vector data(9);
2	fmk	769
272	mhscott	770	data(0) = this->getTag();
		771	data(6) = L;
2	fmk	772
		773	res += theChannel.sendVector(this->getDbTag(), cTag, data);
		774	if (res < 0) {
		775	g3ErrorHandler->warning("%s - failed to send Vector",
		776	"LinearCrdTransf3d::sendSelf");
		777	return res;
		778	}
		779
		780	return res;
		781	}
		782
		783
		784
		785	int
		786	LinearCrdTransf3d::recvSelf(int cTag, Channel &theChannel, FEM_ObjectBroker &theBroker)
		787	{
		788	int res = 0;
		789
272	mhscott	790	static Vector data(9);
2	fmk	791
		792	res += theChannel.recvVector(this->getDbTag(), cTag, data);
		793	if (res < 0) {
		794	g3ErrorHandler->warning("%s - failed to receive Vector",
		795	"LinearCrdTransf3d::recvSelf");
		796	return res;
		797	}
		798
		799	this->setTag((int)data(0));
272	mhscott	800	L = data(6);
2	fmk	801
272	mhscott	802	return res;
2	fmk	803	}
		804
272	mhscott	805
2	fmk	806	const Vector &
		807	LinearCrdTransf3d::getPointGlobalCoordFromLocal(const Vector &xl)
		808	{
		809	static Vector xg(3);
		810
272	mhscott	811	//xg = nodeIPtr->getCrds() + nodeIOffset;
		812	xg = nodeIPtr->getCrds();
		813
		814	if (nodeIOffset) {
		815	xg(0) += nodeIOffset[0];
		816	xg(1) += nodeIOffset[1];
		817	xg(2) += nodeIOffset[2];
		818	}
		819
		820	// xg = xg + Rlj'*xl
		821	//xg.addMatrixTransposeVector(1.0, Rlj, xl, 1.0);
338	mhscott	822	xg(0) += R[0][0]xl(0) + R[1][0]xl(1) + R[2][0]*xl(2);
		823	xg(1) += R[0][1]xl(0) + R[1][1]xl(1) + R[2][1]*xl(2);
		824	xg(2) += R[0][2]xl(0) + R[1][2]xl(1) + R[2][2]*xl(2);
2	fmk	825
		826	return xg;
		827	}
		828
		829
		830	const Vector &
		831	LinearCrdTransf3d::getPointGlobalDisplFromBasic (double xi, const Vector &uxb)
		832	{
		833	// determine global displacements
		834	const Vector &disp1 = nodeIPtr->getTrialDisp();
		835	const Vector &disp2 = nodeJPtr->getTrialDisp();
		836
272	mhscott	837	static double ug[12];
2	fmk	838	for (int i = 0; i < 6; i++)
		839	{
272	mhscott	840	ug[i] = disp1(i);
		841	ug[i+6] = disp2(i);
2	fmk	842	}
		843
		844	// transform global end displacements to local coordinates
272	mhscott	845	//ul.addMatrixVector(0.0, Tlg, ug, 1.0); // ul = Tlg * ug;
		846	static double ul[12];
2	fmk	847
272	mhscott	848	ul[0] = R[0][0]ug[0] + R[0][1]ug[1] + R[0][2]*ug[2];
		849	ul[1] = R[1][0]ug[0] + R[1][1]ug[1] + R[1][2]*ug[2];
		850	ul[2] = R[2][0]ug[0] + R[2][1]ug[1] + R[2][2]*ug[2];
		851
		852	ul[7] = R[1][0]ug[6] + R[1][1]ug[7] + R[1][2]*ug[8];
		853	ul[8] = R[2][0]ug[6] + R[2][1]ug[7] + R[2][2]*ug[8];
		854
		855	static double Wu[3];
		856	if (nodeIOffset) {
		857	Wu[0] = nodeIOffset[2]ug[4] - nodeIOffset[1]ug[5];
		858	Wu[1] = -nodeIOffset[2]ug[3] + nodeIOffset[0]ug[5];
		859	Wu[2] = nodeIOffset[1]ug[3] - nodeIOffset[0]ug[4];
		860
		861	ul[0] += R[0][0]Wu[0] + R[0][1]Wu[1] + R[0][2]*Wu[2];
		862	ul[1] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		863	ul[2] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		864	}
		865
		866	if (nodeJOffset) {
		867	Wu[0] = nodeJOffset[2]ug[10] - nodeJOffset[1]ug[11];
		868	Wu[1] = -nodeJOffset[2]ug[9] + nodeJOffset[0]ug[11];
		869	Wu[2] = nodeJOffset[1]ug[9] - nodeJOffset[0]ug[10];
		870
		871	ul[7] += R[1][0]Wu[0] + R[1][1]Wu[1] + R[1][2]*Wu[2];
		872	ul[8] += R[2][0]Wu[0] + R[2][1]Wu[1] + R[2][2]*Wu[2];
		873	}
		874
2	fmk	875	// compute displacements at point xi, in local coordinates
272	mhscott	876	static double uxl[3];
		877	static Vector uxg(3);
2	fmk	878
272	mhscott	879	uxl[0] = uxb(0) + ul[0];
		880	uxl[1] = uxb(1) + (1-xi)ul[1] + xiul[7];
		881	uxl[2] = uxb(2) + (1-xi)ul[2] + xiul[8];
2	fmk	882
		883	// rotate displacements to global coordinates
272	mhscott	884	// uxg = Rlj'*uxl
		885	//uxg.addMatrixTransposeVector(0.0, Rlj, uxl, 1.0);
		886	uxg(0) = R[0][0]uxl[0] + R[1][0]uxl[1] + R[2][0]*uxl[2];
		887	uxg(1) = R[0][1]uxl[0] + R[1][1]uxl[1] + R[2][1]*uxl[2];
		888	uxg(2) = R[0][2]uxl[0] + R[1][2]uxl[1] + R[2][2]*uxl[2];
2	fmk	889
		890	return uxg;
		891	}
		892
		893
		894	void
		895	LinearCrdTransf3d::Print(ostream &s, int flag)
		896	{
		897	s << "\nCrdTransf: " << this->getTag() << " Type: LinearCrdTransf3d";
272	mhscott	898	if (nodeIOffset)
755	jeremic	899	s << "\tNode I offset: " << nodeIOffset[0] << " " << nodeIOffset[1] << " "<< nodeIOffset[2] << endl;
272	mhscott	900	if (nodeJOffset)
755	jeremic	901	s << "\tNode J offset: " << nodeJOffset[0] << " " << nodeJOffset[1] << " "<< nodeJOffset[2] << endl;
2	fmk	902	}
		903
		904
272	mhscott	905
		906
		907
		908
		909

Subversion Repositories OpenSees

(root)/trunk/SRC/coordTransformation/LinearCrdTransf3d.cpp - Rev 755