Block3D.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: 2004/06/07 23:36:03 $ 00023 // $Source: /usr/local/cvs/OpenSees/SRC/modelbuilder/tcl/Block3D.cpp,v $ 00024 00025 // Written: Ed Love 00026 // Created: 07/01 00027 // 00028 // Description: This file contains the implementation of Block3D. 00029 00030 // 00031 // What: "@(#) Block3D.cpp, revA" 00032 00033 #include <Block3D.h> 00034 00035 00036 //constructor 00037 Block3D::Block3D(int numx, int numy, int numz, 00038 const ID& nodeID, 00039 const Matrix& coorArray ) 00040 : 00041 coor(3), 00042 element(8) 00043 { 00044 this->nx = numx; 00045 this->ny = numy; 00046 this->nz = numz; 00047 00048 this->setUpXl( nodeID, coorArray ); 00049 } 00050 00051 00052 //destructor 00053 Block3D::~Block3D( ) 00054 { } 00055 00056 00057 //set up xl array 00058 void Block3D::setUpXl( const ID &nodeID, const Matrix &coorArray ) 00059 { 00060 00061 int i, j; 00062 00063 for ( i=0; i<8; i++ ){ 00064 if ( nodeID(i) == -1 ) { 00065 opserr << "Warning : in Block3D, block node " 00066 << i 00067 << " is not defined. No Generation will take place." 00068 << endln; 00069 break; 00070 }//end if 00071 }//end for i 00072 00073 00074 //xl = tranpose coorArray(27,3) 00075 for ( i=0; i<3; i++ ) { 00076 for (j=0; j<27; j++ ) 00077 xl[i][j] = coorArray(j,i); 00078 }//end for i 00079 00080 00081 if ( nodeID(8) == -1 ) { 00082 for ( i=0; i<3; i++ ) 00083 xl[i][8] = 0.5*( xl[i][0] + xl[i][4] ); 00084 }//endif 00085 00086 if ( nodeID(9) == -1 ) { 00087 for ( i=0; i<3; i++ ) 00088 xl[i][9] = 0.5*( xl[i][1] + xl[i][5] ); 00089 }//endif 00090 00091 if ( nodeID(10) == -1 ) { 00092 for ( i=0; i<3; i++ ) 00093 xl[i][10] = 0.5*( xl[i][2] + xl[i][6] ); 00094 }//endif 00095 00096 if ( nodeID(11) == -1 ) { 00097 for ( i=0; i<3; i++ ) 00098 xl[i][11] = 0.5*( xl[i][3] + xl[i][7] ); 00099 }//endif 00100 00101 00102 if ( nodeID(12) == -1 ) { 00103 for ( i=0; i<3; i++ ) 00104 xl[i][12] = 0.5*( xl[i][0] + xl[i][1] ); 00105 }//endif 00106 00107 if ( nodeID(13) == -1 ) { 00108 for ( i=0; i<3; i++ ) 00109 xl[i][13] = 0.5*( xl[i][1] + xl[i][2] ); 00110 }//endif 00111 00112 if ( nodeID(14) == -1 ) { 00113 for ( i=0; i<3; i++ ) 00114 xl[i][14] = 0.5*( xl[i][2] + xl[i][3] ); 00115 }//endif 00116 00117 if ( nodeID(15) == -1 ) { 00118 for ( i=0; i<3; i++ ) 00119 xl[i][15] = 0.5*( xl[i][0] + xl[i][3] ); 00120 }//endif 00121 00122 00123 if ( nodeID(16) == -1 ) { 00124 for ( i=0; i<3; i++ ) 00125 xl[i][16] = 0.25*( xl[i][0] + xl[i][1] + xl[i][2] + xl[i][3] ); 00126 }//endif 00127 00128 00129 if ( nodeID(17) == -1 ) { 00130 for ( i=0; i<3; i++ ) 00131 xl[i][17] = 0.5*( xl[i][4] + xl[i][5] ); 00132 }//endif 00133 00134 if ( nodeID(18) == -1 ) { 00135 for ( i=0; i<3; i++ ) 00136 xl[i][18] = 0.5*( xl[i][5] + xl[i][6] ); 00137 }//endif 00138 00139 if ( nodeID(19) == -1 ) { 00140 for ( i=0; i<3; i++ ) 00141 xl[i][19] = 0.5*( xl[i][6] + xl[i][7] ); 00142 }//endif 00143 00144 if ( nodeID(20) == -1 ) { 00145 for ( i=0; i<3; i++ ) 00146 xl[i][20] = 0.5*( xl[i][4] + xl[i][7] ); 00147 }//endif 00148 00149 00150 if ( nodeID(21) == -1 ) { 00151 for ( i=0; i<3; i++ ) 00152 xl[i][21] = 0.25*( xl[i][4] + xl[i][5] + xl[i][6] + xl[i][7] ); 00153 }//endif 00154 00155 00156 if ( nodeID(22) == -1 ) { 00157 for ( i=0; i<3; i++ ) 00158 xl[i][22] = 0.25*( xl[i][0] + xl[i][1] + xl[i][5] + xl[i][4] ); 00159 }//endif 00160 00161 if ( nodeID(23) == -1 ) { 00162 for ( i=0; i<3; i++ ) 00163 xl[i][23] = 0.25*( xl[i][1] + xl[i][2] + xl[i][6] + xl[i][5] ); 00164 }//endif 00165 00166 00167 if ( nodeID(24) == -1 ) { 00168 for ( i=0; i<3; i++ ) 00169 xl[i][24] = 0.25*( xl[i][3] + xl[i][2] + xl[i][6] + xl[i][7] ); 00170 }//endif 00171 00172 if ( nodeID(25) == -1 ) { 00173 for ( i=0; i<3; i++ ) 00174 xl[i][25] = 0.25*( xl[i][0] + xl[i][3] + xl[i][7] + xl[i][4] ); 00175 }//endif 00176 00177 00178 00179 if ( nodeID(26) == -1 ) { 00180 for ( i=0; i<3; i++ ) 00181 xl[i][26] = 0.125*( xl[i][0] + xl[i][1] + xl[i][2] + xl[i][3] + 00182 xl[i][4] + xl[i][5] + xl[i][6] + xl[i][7] ); 00183 }//endif 00184 00185 return; 00186 } 00187 00188 00189 //generate node 00190 const Vector& 00191 Block3D::getNodalCoords( int i, int j, int k ) 00192 { 00193 00194 /* loop as follows (in pseudocode) 00195 for ( k = 0, nz ) { 00196 for ( j = 0, ny ) { 00197 for ( i = 0, nx ) 00198 call getNodalCoords(i,j,k); 00199 } 00200 } 00201 */ 00202 00203 double hx = 2.0 / nx; 00204 00205 double hy = 2.0 / ny; 00206 00207 double hhz = 2.0 / nz; 00208 00209 double x = -1.0 + (i*hx); 00210 00211 double y = -1.0 + (j*hy); 00212 00213 double z = -1.0 + (k*hhz); 00214 00215 coor(0) = x; 00216 coor(1) = y; 00217 coor(2) = z; 00218 00219 this->transformNodalCoordinates( ); 00220 00221 return coor; 00222 } 00223 00224 00225 //generate element 00226 const ID& 00227 Block3D::getElementNodes( int i, int j, int k ) 00228 { 00229 00230 int nenx = nx + 1; 00231 00232 int neny = ny + 1; 00233 00234 int nInXYplane = nenx * neny; 00235 00236 00237 int node1, node2, node3, node4; 00238 int node5, node6, node7, node8; 00239 00240 00241 node1 = i + (j*nenx) + (k*nInXYplane); 00242 node2 = node1 + 1; 00243 node3 = node2 + nenx; 00244 node4 = node1 + nenx; 00245 00246 node5 = node1 + nInXYplane; 00247 node6 = node2 + nInXYplane; 00248 node7 = node3 + nInXYplane; 00249 node8 = node4 + nInXYplane; 00250 00251 element(0) = node1; 00252 element(1) = node2; 00253 element(2) = node3; 00254 element(3) = node4; 00255 00256 element(4) = node5; 00257 element(5) = node6; 00258 element(6) = node7; 00259 element(7) = node8; 00260 00261 return element; 00262 } 00263 00264 00265 00266 //transform to real coordinates 00267 void Block3D::transformNodalCoordinates( ) 00268 { 00269 00270 static double shape[27]; 00271 00272 static double natCoor[3]; 00273 00274 int j, dim; 00275 00276 natCoor[0] = coor(0); 00277 natCoor[1] = coor(1); 00278 natCoor[2] = coor(2); 00279 00280 coor.Zero( ); 00281 00282 this->shape3d( natCoor[0], natCoor[1], natCoor[2], shape ); 00283 00284 for ( j=0; j<27; j++ ) { 00285 00286 for ( dim=0; dim<3; dim++ ) 00287 coor(dim) += shape[j]*xl[dim][j]; 00288 00289 } //end for j 00290 00291 return; 00292 00293 } 00294 00295 00296 00297 //shape functions 00298 void Block3D::shape3d( double r, double s, double t, 00299 double shape[27] ) 00300 { 00301 00302 static const int ri[] = {-1, 1, 1,-1, -1, 1, 1,-1, -1, 1, 1,-1, 0, 1, 0,-1, 0, 0, 1, 0,-1, 0, 0, 1, 0,-1, 0}; 00303 00304 static const int si[] = {-1,-1, 1, 1, -1,-1, 1, 1, -1,-1, 1, 1, -1, 0, 1, 0, 0, -1, 0, 1, 0, 0, -1, 0, 1, 0, 0}; 00305 00306 static const int ti[] = {-1,-1,-1,-1, 1, 1, 1, 1, 0, 0, 0, 0, -1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0}; 00307 00308 00309 static const double d1 = 1.0; 00310 static const double d2 = 0.5; // = 1.0/2.0; 00311 static const double d4 = 0.25; // = 1.0/4.0; 00312 static const double d8 = 0.125; // = 1.0/8.0; 00313 00314 00315 double rr = r*r; 00316 double ss = s*s; 00317 double tt = t*t; 00318 00319 double r0, s0, t0; 00320 int k, kk; 00321 00322 //shape functions for 27-node element 00323 for ( k=1; k<=27; k++ ) { 00324 00325 kk = k-1; //C-style numbering 00326 00327 r0 = r*ri[kk]; 00328 s0 = s*si[kk]; 00329 t0 = t*ti[kk]; 00330 00331 //corner nodes top/bottom 00332 if ( k>=1 && k<=8 ) 00333 shape[kk] = d8*(rr+r0) *(ss+s0) *(tt+t0); 00334 00335 //corner nodes midside 00336 if ( k>=9 && k<=12 ) 00337 shape[kk] = d4*(rr+r0) *(ss+s0) *(d1-tt); 00338 00339 //midside nodes top/bottom r-dir 00340 if ( k==13 || k==15 || k==18 || k==20 ) 00341 shape[kk] = d4*(d1-rr)*(ss+s0) *(tt+t0); 00342 00343 //midside nodes top/bottom s-dir 00344 if ( k==14 || k==16 || k==19 || k==21 ) 00345 shape[kk] = d4*(rr+r0) *(d1-ss)*(tt+t0); 00346 00347 //midside nodes mid plane r-dir 00348 if ( k==23 || k==25 ) 00349 shape[kk] = d2*(d1-rr)*(ss+s0) *(d1-tt); 00350 00351 //midside nodes mid plane s-dir 00352 if ( k==24 || k==26 ) 00353 shape[kk] = d2*(rr+r0) *(d1-ss)*(d1-tt); 00354 00355 //central nodes top/bottom 00356 if ( k==17 || k==22 ) 00357 shape[kk] = d2*(d1-rr)*(d1-ss)*(tt+t0); 00358 00359 if ( k==27 ) 00360 shape[kk] = (d1-rr)*(d1-ss)*(d1-tt); 00361 00362 }//end for k 00363 00364 return; 00365 } 00366 00367 00368 /* 00369 subroutine shp04(shp,glu,glo,gu,eu,to,xjac,detj,r,s,t,xl,ul) 00370 c----------------------------------------------------------------------- 00371 c.....compute shape functions and their derivatives for linear,quadratic 00372 c.....lagrangian and serendipity isoparametric 3-d elements 00373 c.....global coordinate system x,y,z 00374 c.....local coordinate system xsi,eta,zeta 00375 c----------------------------------------------------------------------- 00376 implicit double precision (a-h,o-z) 00377 dimension shp(4,27),ri(27),si(27),ti(27) 00378 dimension xt(3,3),cro(3),to(3,3), 00379 + xl(3,27),ul(3,27),xu(3,27),xjac(3,3), 00380 + glu(3,3),glo(3,3),gu(3,3),gko(3,3),gku(3,3),eu(3,3) 00381 common /eldata/ dm,n,ma,mct,iel,nel 00382 common /iofile/ ior,iow 00383 data ri /-1, 1, 1,-1, -1, 1, 1,-1, -1, 1, 1,-1, 0, 1, 0,-1, 0, 00384 1 0, 1, 0,-1, 0, 0, 1, 0,-1, 0/ 00385 data si /-1,-1, 1, 1, -1,-1, 1, 1, -1,-1, 1, 1, -1, 0, 1, 0, 0, 00386 1 -1, 0, 1, 0, 0, -1, 0, 1, 0, 0/ 00387 data ti /-1,-1,-1,-1, 1, 1, 1, 1, 0, 0, 0, 0, -1,-1,-1,-1,-1, 00388 1 1, 1, 1, 1, 1, 0, 0, 0, 0, 0/ 00389 00390 d1=1.0d0 00391 d2=1.0d0/2.0d0 00392 d8=1.0d0/8.0d0 00393 d4=1.0d0/4.0d0 00394 dz=2.0d0 00395 rr=r*r 00396 ss=s*s 00397 tt=t*t 00398 00399 00400 c.... shape functions for 27-node element 00401 do l = 1,27 00402 r0 = r*ri(l) 00403 s0 = s*si(l) 00404 t0 = t*ti(l) 00405 c corner nodes top/bottom 00406 if(l.ge.1.and.l.le.8) then 00407 shp(4,l) = d8*(rr+r0) *(ss+s0) *(tt+t0) 00408 c corner nodes mid-side 00409 elseif(l.ge.9.and.l.le.12) then 00410 shp(4,l) = d4*(rr+r0) *(ss+s0) *(d1-tt) 00411 c midside nodes top/bottom r-dir 00412 elseif(l.eq.13.or.l.eq.15.or.l.eq.18.or.l.eq.20) then 00413 shp(4,l) = d4*(d1-rr)*(ss+s0) *(tt+t0) 00414 c midside nodes top/bottom s-dir 00415 elseif(l.eq.14.or.l.eq.16.or.l.eq.19.or.l.eq.21) then 00416 shp(4,l) = d4*(rr+r0) *(d1-ss)*(tt+t0) 00417 c midside nodes mid plane r-dir 00418 elseif(l.eq.23.or.l.eq.25) then 00419 shp(4,l) = d2*(d1-rr)*(ss+s0) *(d1-tt) 00420 c midside nodes mid plane s-dir 00421 elseif(l.eq.24.or.l.eq.26) then 00422 shp(4,l) = d2*(rr+r0) *(d1-ss)*(d1-tt) 00423 c central nodes top/bottom 00424 elseif(l.eq.17.or.l.eq.22) then 00425 shp(4,l) = d2*(d1-rr)*(d1-ss)*(tt+t0) 00426 elseif(l.eq.27) then 00427 shp(4,l) = (d1-rr)*(d1-ss)*(d1-tt) 00428 endif 00429 enddo 00430 */
Generated on Mon Oct 23 15:05:23 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents