R3vectors.cppGo to the documentation of this file.00001 /* ****************************************************************** ** 00002 ** OpenSees - Open System for Earthquake Engineering Simulation ** 00003 ** Pacific Earthquake Engineering Research Center ** 00004 ** ** 00005 ** ** 00006 ** (C) Copyright 1999, The Regents of the University of California ** 00007 ** All Rights Reserved. ** 00008 ** ** 00009 ** Commercial use of this program without express permission of the ** 00010 ** University of California, Berkeley, is strictly prohibited. See ** 00011 ** file 'COPYRIGHT' in main directory for information on usage and ** 00012 ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** 00013 ** ** 00014 ** Developed by: ** 00015 ** Frank McKenna (fmckenna@ce.berkeley.edu) ** 00016 ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** 00017 ** Filip C. Filippou (filippou@ce.berkeley.edu) ** 00018 ** ** 00019 ** ****************************************************************** */ 00020 00021 // $Revision: 1.3 $ 00022 // $Date: 2002/12/05 22:20:44 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/element/shell/R3vectors.cpp,v $ 00024 00025 // Ed "C++" Love 00026 00027 #include <R3vectors.h> 00028 #include <Vector.h> 00029 #include <Matrix.h> 00030 #include <math.h> 00031 00032 #define sign(a) ( (a)>0 ? 1:-1 ) 00033 00034 double LovelyInnerProduct( const Vector &v, const Vector &w ) 00035 { 00036 int i ; 00037 00038 double dot = 0.0 ; 00039 00040 for ( i = 0; i < 3; i++ ) 00041 dot += v(i)*w(i) ; 00042 00043 return dot ; 00044 00045 } 00046 00047 00048 double LovelyNorm( const Vector &v ) 00049 { 00050 return sqrt( LovelyInnerProduct(v,v) ) ; 00051 } 00052 00053 00054 Vector LovelyCrossProduct( const Vector &v, const Vector &w ) 00055 { 00056 00057 Vector cross(3) ; 00058 00059 cross(0) = v(1)*w(2) - v(2)*w(1) ; 00060 00061 cross(1) = v(2)*w(0) - v(0)*w(2) ; 00062 00063 cross(2) = v(0)*w(1) - v(1)*w(0) ; 00064 00065 return cross ; 00066 00067 } 00068 00069 00070 Vector LovelyEig( const Matrix &M ) 00071 { 00072 //.... compute eigenvalues and vectors for a 3 x 3 symmetric matrix 00073 // 00074 //.... INPUTS: 00075 // M(3,3) - matrix with initial values (only upper half used) 00076 // 00077 //.... OUTPUTS 00078 // v(3,3) - matrix of eigenvectors (by column) 00079 // d(3) - eigenvalues associated with columns of v 00080 // rot - number of rotations to diagonalize 00081 // 00082 //---------------------------------------------------------------eig3== 00083 00084 //.... Storage done as follows: 00085 // 00086 // | v(1,1) v(1,2) v(1,3) | | d(1) a(1) a(3) | 00087 // | v(2,1) v(2,2) v(2,3) | = | a(1) d(2) a(2) | 00088 // | v(3,1) v(3,2) v(3,3) | | a(3) a(2) d(3) | 00089 // 00090 // Transformations performed on d(i) and a(i) and v(i,j) become 00091 // the eigenvectors. 00092 // 00093 //---------------------------------------------------------------eig3== 00094 00095 int rot, its, i, j , k ; 00096 double g, h, aij, sm, thresh, t, c, s, tau ; 00097 00098 static Matrix v(3,3) ; 00099 static Vector d(3) ; 00100 static Vector a(3) ; 00101 static Vector b(3) ; 00102 static Vector z(3) ; 00103 00104 static const double tol = 1.0e-08 ; 00105 00106 // set v = M 00107 v = M ; 00108 00109 //.... move array into one-d arrays 00110 00111 a(0) = v(0,1) ; 00112 a(1) = v(1,2) ; 00113 a(2) = v(2,0) ; 00114 00115 00116 for ( i = 0; i < 3; i++ ) { 00117 d(i) = v(i,i) ; 00118 b(i) = v(i,i) ; 00119 z(i) = 0.0 ; 00120 00121 for ( j = 0; j < 3; j++ ) 00122 v(i,j) = 0.0 ; 00123 00124 v(i,i) = 1.0 ; 00125 00126 } //end for i 00127 00128 rot = 0 ; 00129 its = 0 ; 00130 00131 sm = fabs(a(0)) + fabs(a(1)) + fabs(a(2)) ; 00132 00133 while ( sm > tol ) { 00134 //.... set convergence test and threshold 00135 if ( its < 3 ) 00136 thresh = 0.011*sm ; 00137 else 00138 thresh = 0.0 ; 00139 00140 //.... perform sweeps for rotations 00141 for ( i = 0; i < 3; i++ ) { 00142 00143 j = (i+1)%3; 00144 k = (j+1)%3; 00145 00146 aij = a(i) ; 00147 00148 g = 100.0 * fabs(aij) ; 00149 00150 if ( fabs(d(i)) + g != fabs(d(i)) || 00151 fabs(d(j)) + g != fabs(d(j)) ) { 00152 00153 if ( fabs(aij) > thresh ) { 00154 00155 a(i) = 0.0 ; 00156 h = d(j) - d(i) ; 00157 00158 if( fabs(h)+g == fabs(h) ) 00159 t = aij / h ; 00160 else { 00161 //t = 2.0 * sign(h/aij) / ( fabs(h/aij) + sqrt(4.0+(h*h/aij/aij))); 00162 double hDIVaij = h/aij; 00163 if (hDIVaij > 0.0) 00164 t = 2.0 / ( hDIVaij + sqrt(4.0+(hDIVaij * hDIVaij))); 00165 else 00166 t = - 2.0 / (-hDIVaij + sqrt(4.0+(hDIVaij * hDIVaij))); 00167 } 00168 00169 //.... set rotation parameters 00170 00171 c = 1.0 / sqrt(1.0 + t*t) ; 00172 s = t*c ; 00173 tau = s / (1.0 + c) ; 00174 00175 //.... rotate diagonal terms 00176 00177 h = t * aij ; 00178 z(i) = z(i) - h ; 00179 z(j) = z(j) + h ; 00180 d(i) = d(i) - h ; 00181 d(j) = d(j) + h ; 00182 00183 //.... rotate off-diagonal terms 00184 00185 h = a(j) ; 00186 g = a[k] ; 00187 a(j) = h + s*(g - h*tau) ; 00188 a(k) = g - s*(h + g*tau) ; 00189 00190 //.... rotate eigenvectors 00191 00192 for ( k = 0; k < 3; k++ ) { 00193 g = v(k,i) ; 00194 h = v(k,j) ; 00195 v(k,i) = g - s*(h + g*tau) ; 00196 v(k,j) = h + s*(g - h*tau) ; 00197 } // end for k 00198 00199 rot = rot + 1 ; 00200 00201 } // end if fabs > thresh 00202 } //else 00203 else 00204 a(i) = 0.0 ; 00205 00206 } // end for i 00207 00208 //.... update the diagonal terms 00209 for ( i = 0; i < 3; i++ ) { 00210 b(i) = b(i) + z(i) ; 00211 d(i) = b(i) ; 00212 z(i) = 0.0 ; 00213 } // end for i 00214 00215 its += 1 ; 00216 00217 sm = fabs(a(0)) + fabs(a(1)) + fabs(a(2)) ; 00218 00219 } //end while sm 00220 00221 return d ; 00222 } 00223
Generated on Mon Oct 23 15:05:09 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents