OpenSees » SRC » element » brick

shp3d.cpp

Go to the documentation of this file.
00001 
00002 // $Revision: 1.1 $
00003 // $Date: 2001/07/11 21:54:41 $
00004 // $Source: /usr/local/cvs/OpenSees/SRC/element/brick/shp3d.cpp,v $
00005 
00006 // Ed "C++" Love
00007 //
00008 // 3-d isoparametric 8-node element shape function
00009 //
00010 
00011 /* 
00012       Purpose: Compute 3-d isoparametric 8-node element shape
00013                functions and their derivatives w/r x,y,z
00014 
00015       Inputs:
00016          ss[3]     - Natural coordinates of point
00017          xl[3][8]  - Nodal coordinates for element
00018 
00019       Outputs:
00020          xsj        - Jacobian determinant at point
00021          shp[4][8]  - Shape functions and derivatives at point
00022                      shp[0][i] = dN_i/dx
00023                      shp[1][i] = dN_i/dy
00024                      shp[2][i] = dN_i/dzc
00025                      shp[3][i] =  N_i
00026 */
00027 
00028 void  shp3d( const double ss[3],
00029              double &xsj,
00030              double shp[4][8],
00031              const double xl[3][8]   )
00032 {
00033     int   i, j, k ;
00034 
00035     double rxsj, ap1, am1, ap2, am2, ap3, am3, c1,c2,c3 ;
00036 
00037     static double xs[3][3] ; 
00038     static double ad[3][3] ;
00039 
00040 
00041       //Compute shape functions and their natural coord. derivatives
00042 
00043       ap1 = 1.0 + ss[0] ;
00044       am1 = 1.0 - ss[0] ;
00045       ap2 = 1.0 + ss[1] ;
00046       am2 = 1.0 - ss[1] ;
00047       ap3 = 1.0 + ss[2] ;
00048       am3 = 1.0 - ss[2] ;
00049 
00050       //Compute for ( - , - ) values
00051 
00052       c1      = 0.125*am1*am2 ;
00053       c2      = 0.125*am2*am3 ;
00054       c3      = 0.125*am1*am3 ;
00055       shp[0][0] = -c2 ;
00056       shp[0][1] =  c2 ;
00057       shp[1][0] = -c3 ;
00058       shp[1][3] =  c3 ;
00059       shp[2][0] = -c1 ;
00060       shp[2][4] =  c1 ;
00061       shp[3][0] =  c1*am3 ;
00062       shp[3][4] =  c1*ap3 ;
00063 
00064       //Compute for ( + , + ) values
00065 
00066       c1      = 0.125*ap1*ap2 ;
00067       c2      = 0.125*ap2*ap3 ;
00068       c3      = 0.125*ap1*ap3 ;
00069       shp[0][7] = -c2 ;
00070       shp[0][6] =  c2 ;
00071       shp[1][5] = -c3 ;
00072       shp[1][6] =  c3 ;
00073       shp[2][2] = -c1 ;
00074       shp[2][6] =  c1 ;
00075       shp[3][2] =  c1*am3 ;
00076       shp[3][6] =  c1*ap3 ;
00077 
00078       //Compute for ( - , + ) values
00079 
00080       c1      = 0.125*am1*ap2 ;
00081       c2      = 0.125*am2*ap3 ;
00082       c3      = 0.125*am1*ap3 ;
00083       shp[0][4] = -c2 ;
00084       shp[0][5] =  c2 ; 
00085       shp[1][4] = -c3 ;
00086       shp[1][7] =  c3 ;
00087       shp[2][3] = -c1 ;
00088       shp[2][7] =  c1 ;
00089       shp[3][3] =  c1*am3 ;
00090       shp[3][7] =  c1*ap3 ;
00091 
00092       //Compute for ( + , - ) values
00093 
00094       c1      = 0.125*ap1*am2 ;
00095       c2      = 0.125*ap2*am3 ;
00096       c3      = 0.125*ap1*am3 ;
00097       shp[0][3] = -c2 ;
00098       shp[0][2] =  c2 ;
00099       shp[1][1] = -c3 ;
00100       shp[1][2] =  c3 ;
00101       shp[2][1] = -c1 ;
00102       shp[2][5] =  c1 ;
00103       shp[3][1] =  c1*am3 ;
00104       shp[3][5] =  c1*ap3 ;
00105 
00106       //Compute jacobian transformation
00107 
00108       for ( j=0; j<3; j++ ) {
00109 
00110         xs[j][0] = ( xl[j][1] - xl[j][0] )*shp[0][1]
00111                  + ( xl[j][2] - xl[j][3] )*shp[0][2]
00112                  + ( xl[j][5] - xl[j][4] )*shp[0][5]
00113                  + ( xl[j][6] - xl[j][7] )*shp[0][6] ;
00114 
00115         xs[j][1] = ( xl[j][2] - xl[j][1] )*shp[1][2]
00116                  + ( xl[j][3] - xl[j][0] )*shp[1][3]
00117                  + ( xl[j][6] - xl[j][5] )*shp[1][6]
00118                  + ( xl[j][7] - xl[j][4] )*shp[1][7] ;
00119 
00120         xs[j][2] = ( xl[j][4] - xl[j][0] )*shp[2][4]
00121                  + ( xl[j][5] - xl[j][1] )*shp[2][5]
00122                  + ( xl[j][6] - xl[j][2] )*shp[2][6]
00123                  + ( xl[j][7] - xl[j][3] )*shp[2][7] ;
00124 
00125       } //end for j     
00126 
00127 
00128       //Compute adjoint to jacobian
00129 
00130       ad[0][0] = xs[1][1]*xs[2][2] - xs[1][2]*xs[2][1] ;
00131       ad[0][1] = xs[2][1]*xs[0][2] - xs[2][2]*xs[0][1] ;
00132       ad[0][2] = xs[0][1]*xs[1][2] - xs[0][2]*xs[1][1] ;
00133 
00134       ad[1][0] = xs[1][2]*xs[2][0] - xs[1][0]*xs[2][2] ;
00135       ad[1][1] = xs[2][2]*xs[0][0] - xs[2][0]*xs[0][2] ;
00136       ad[1][2] = xs[0][2]*xs[1][0] - xs[0][0]*xs[1][2] ;
00137 
00138       ad[2][0] = xs[1][0]*xs[2][1] - xs[1][1]*xs[2][0] ;
00139       ad[2][1] = xs[2][0]*xs[0][1] - xs[2][1]*xs[0][0] ; 
00140       ad[2][2] = xs[0][0]*xs[1][1] - xs[0][1]*xs[1][0] ;
00141 
00142       //Compute determinant of jacobian
00143 
00144       xsj  = xs[0][0]*ad[0][0] + xs[0][1]*ad[1][0] + xs[0][2]*ad[2][0] ;
00145       rxsj = 1.0/xsj ;
00146 
00147       //Compute jacobian inverse
00148 
00149       for ( j=0; j<3; j++ ) {
00150 
00151         for ( i=0; i<3; i++ ) 
00152              xs[i][j] = ad[i][j]*rxsj ;
00153 
00154       } //end for j
00155 
00156 
00157       //Compute derivatives with repect to global coords.
00158 
00159       for ( k=0; k<8; k++) {
00160 
00161         c1 = shp[0][k]*xs[0][0] + shp[1][k]*xs[1][0] + shp[2][k]*xs[2][0] ;
00162         c2 = shp[0][k]*xs[0][1] + shp[1][k]*xs[1][1] + shp[2][k]*xs[2][1] ;
00163         c3 = shp[0][k]*xs[0][2] + shp[1][k]*xs[1][2] + shp[2][k]*xs[2][2] ;
00164 
00165         shp[0][k] = c1 ;
00166         shp[1][k] = c2 ;
00167         shp[2][k] = c3 ;
00168 
00169       } //end for k
00170 
00171    return ;
00172 }

Generated on Mon Oct 23 15:05:05 2006 for OpenSees by doxygen 1.5.0