MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /*========================================================================= 00002 00003 Module: $RCSfile: V_GaussIntegration.cpp,v $ 00004 00005 Copyright (c) 2006 Sandia Corporation. 00006 All rights reserved. 00007 See Copyright.txt or http://www.kitware.com/Copyright.htm for details. 00008 00009 This software is distributed WITHOUT ANY WARRANTY; without even 00010 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00011 PURPOSE. See the above copyright notice for more information. 00012 00013 =========================================================================*/ 00014 00015 /* 00016 * 00017 * GaussIntegration.cpp performs Gauss Integrations 00018 * 00019 * This file is part of VERDICT 00020 * 00021 */ 00022 00023 #define VERDICT_EXPORTS 00024 00025 #include "moab/verdict.h" 00026 #include "V_GaussIntegration.hpp" 00027 00028 #include <math.h> 00029 00030 double verdictSqrt2 = sqrt( (double)2.0 ); 00031 00032 int numberGaussPoints; 00033 int numberNodes; 00034 int numberDims; 00035 double gaussPointY[maxNumberGaussPoints]; 00036 double gaussWeight[maxNumberGaussPoints]; 00037 double shapeFunction[maxTotalNumberGaussPoints][maxNumberNodes]; 00038 double dndy1GaussPts[maxTotalNumberGaussPoints][maxNumberNodes]; 00039 double dndy2GaussPts[maxTotalNumberGaussPoints][maxNumberNodes]; 00040 double dndy3GaussPts[maxTotalNumberGaussPoints][maxNumberNodes]; 00041 double totalGaussWeight[maxTotalNumberGaussPoints]; 00042 int totalNumberGaussPts; 00043 double y1Area[maxNumberGaussPointsTri]; 00044 double y2Area[maxNumberGaussPointsTri]; 00045 double y1Volume[maxNumberGaussPointsTet]; 00046 double y2Volume[maxNumberGaussPointsTet]; 00047 double y3Volume[maxNumberGaussPointsTet]; 00048 double y4Volume[maxNumberGaussPointsTet]; 00049 00050 void GaussIntegration::initialize( int n, int m, int dim, int tri ) 00051 { 00052 numberGaussPoints = n; 00053 numberNodes = m; 00054 numberDims = dim; 00055 00056 if( tri == 1 ) 00057 // triangular element 00058 { 00059 if( numberDims == 2 ) 00060 totalNumberGaussPts = numberGaussPoints; 00061 else if( numberDims == 3 ) 00062 totalNumberGaussPts = numberGaussPoints; 00063 } 00064 else if( tri == 0 ) 00065 { 00066 if( numberDims == 2 ) 00067 totalNumberGaussPts = numberGaussPoints * numberGaussPoints; 00068 else if( numberDims == 3 ) 00069 totalNumberGaussPts = numberGaussPoints * numberGaussPoints * numberGaussPoints; 00070 } 00071 } 00072 00073 void GaussIntegration::get_shape_func( double shape_function[], double dndy1_at_gauss_pts[], 00074 double dndy2_at_gauss_pts[], double gauss_weight[] ) 00075 { 00076 int i, j; 00077 for( i = 0; i < totalNumberGaussPts; i++ ) 00078 { 00079 for( j = 0; j < numberNodes; j++ ) 00080 { 00081 shape_function[i * maxNumberNodes + j] = shapeFunction[i][j]; 00082 dndy1_at_gauss_pts[i * maxNumberNodes + j] = dndy1GaussPts[i][j]; 00083 dndy2_at_gauss_pts[i * maxNumberNodes + j] = dndy2GaussPts[i][j]; 00084 } 00085 } 00086 00087 for( i = 0; i < totalNumberGaussPts; i++ ) 00088 gauss_weight[i] = totalGaussWeight[i]; 00089 } 00090 00091 void GaussIntegration::get_shape_func( double shape_function[], double dndy1_at_gauss_pts[], 00092 double dndy2_at_gauss_pts[], double dndy3_at_gauss_pts[], double gauss_weight[] ) 00093 { 00094 int i, j; 00095 for( i = 0; i < totalNumberGaussPts; i++ ) 00096 { 00097 for( j = 0; j < numberNodes; j++ ) 00098 { 00099 shape_function[i * maxNumberNodes + j] = shapeFunction[i][j]; 00100 dndy1_at_gauss_pts[i * maxNumberNodes + j] = dndy1GaussPts[i][j]; 00101 dndy2_at_gauss_pts[i * maxNumberNodes + j] = dndy2GaussPts[i][j]; 00102 dndy3_at_gauss_pts[i * maxNumberNodes + j] = dndy3GaussPts[i][j]; 00103 } 00104 } 00105 00106 for( i = 0; i < totalNumberGaussPts; i++ ) 00107 gauss_weight[i] = totalGaussWeight[i]; 00108 } 00109 00110 void GaussIntegration::get_gauss_pts_and_weight() 00111 { 00112 00113 switch( numberGaussPoints ) 00114 { 00115 case 1: 00116 gaussPointY[0] = 0.0; 00117 gaussWeight[0] = 2.0; 00118 break; 00119 case 2: 00120 gaussPointY[0] = -0.577350269189626; 00121 gaussPointY[1] = 0.577350269189626; 00122 gaussWeight[0] = 1.0; 00123 gaussWeight[1] = 1.0; 00124 break; 00125 case 3: 00126 gaussPointY[0] = -0.774596669241483; 00127 gaussPointY[1] = 0.0; 00128 gaussPointY[2] = 0.774596669241483; 00129 gaussWeight[0] = 0.555555555555555; 00130 gaussWeight[1] = 0.888888888888889; 00131 gaussWeight[2] = 0.555555555555555; 00132 break; 00133 } 00134 } 00135 00136 void GaussIntegration::calculate_shape_function_2d_quad() 00137 { 00138 int ife = 0, i, j; 00139 double y1, y2; 00140 get_gauss_pts_and_weight(); 00141 00142 switch( numberNodes ) 00143 { 00144 case 4: 00145 for( i = 0; i < numberGaussPoints; i++ ) 00146 { 00147 for( j = 0; j < numberGaussPoints; j++ ) 00148 { 00149 y1 = gaussPointY[i]; 00150 y2 = gaussPointY[j]; 00151 shapeFunction[ife][0] = 0.25 * ( 1 - y1 ) * ( 1 - y2 ); 00152 shapeFunction[ife][1] = 0.25 * ( 1 + y1 ) * ( 1 - y2 ); 00153 shapeFunction[ife][2] = 0.25 * ( 1 + y1 ) * ( 1 + y2 ); 00154 shapeFunction[ife][3] = 0.25 * ( 1 - y1 ) * ( 1 + y2 ); 00155 00156 dndy1GaussPts[ife][0] = -0.25 * ( 1 - y2 ); 00157 dndy1GaussPts[ife][1] = 0.25 * ( 1 - y2 ); 00158 dndy1GaussPts[ife][2] = 0.25 * ( 1 + y2 ); 00159 dndy1GaussPts[ife][3] = -0.25 * ( 1 + y2 ); 00160 00161 dndy2GaussPts[ife][0] = -0.25 * ( 1 - y1 ); 00162 dndy2GaussPts[ife][1] = -0.25 * ( 1 + y1 ); 00163 dndy2GaussPts[ife][2] = 0.25 * ( 1 + y1 ); 00164 dndy2GaussPts[ife][3] = 0.25 * ( 1 - y1 ); 00165 00166 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j]; 00167 ife++; 00168 } 00169 } 00170 break; 00171 case 8: 00172 for( i = 0; i < numberGaussPoints; i++ ) 00173 { 00174 for( j = 0; j < numberGaussPoints; j++ ) 00175 { 00176 y1 = gaussPointY[i]; 00177 y2 = gaussPointY[j]; 00178 shapeFunction[ife][0] = 0.25 * ( 1 - y1 ) * ( 1 - y2 ) * ( -y1 - y2 - 1 ); 00179 shapeFunction[ife][1] = 0.25 * ( 1 + y1 ) * ( 1 - y2 ) * ( y1 - y2 - 1 ); 00180 shapeFunction[ife][2] = 0.25 * ( 1 + y1 ) * ( 1 + y2 ) * ( y1 + y2 - 1 ); 00181 shapeFunction[ife][3] = 0.25 * ( 1 - y1 ) * ( 1 + y2 ) * ( -y1 + y2 - 1 ); 00182 shapeFunction[ife][4] = 0.5 * ( 1 - y1 * y1 ) * ( 1 - y2 ); 00183 shapeFunction[ife][5] = 0.5 * ( 1 - y2 * y2 ) * ( 1 + y1 ); 00184 shapeFunction[ife][6] = 0.5 * ( 1 - y1 * y1 ) * ( 1 + y2 ); 00185 shapeFunction[ife][7] = 0.5 * ( 1 - y2 * y2 ) * ( 1 - y1 ); 00186 00187 dndy1GaussPts[ife][0] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 + y2 ); 00188 dndy1GaussPts[ife][1] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 - y2 ); 00189 dndy1GaussPts[ife][2] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 + y2 ); 00190 dndy1GaussPts[ife][3] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 - y2 ); 00191 00192 dndy1GaussPts[ife][4] = -y1 * ( 1 - y2 ); 00193 dndy1GaussPts[ife][5] = 0.5 * ( 1 - y2 * y2 ); 00194 dndy1GaussPts[ife][6] = -y1 * ( 1 + y2 ); 00195 dndy1GaussPts[ife][7] = -0.5 * ( 1 - y2 * y2 ); 00196 00197 dndy2GaussPts[ife][0] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 + y1 ); 00198 dndy2GaussPts[ife][1] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 - y1 ); 00199 dndy2GaussPts[ife][2] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 + y1 ); 00200 dndy2GaussPts[ife][3] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 - y1 ); 00201 00202 dndy2GaussPts[ife][4] = -0.5 * ( 1 - y1 * y1 ); 00203 dndy2GaussPts[ife][5] = -y2 * ( 1 + y1 ); 00204 dndy2GaussPts[ife][6] = 0.5 * ( 1 - y1 * y1 ); 00205 dndy2GaussPts[ife][7] = -y2 * ( 1 - y1 ); 00206 00207 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j]; 00208 ife++; 00209 } 00210 } 00211 break; 00212 } 00213 } 00214 00215 void GaussIntegration::calculate_shape_function_3d_hex() 00216 { 00217 int ife = 0, i, j, k, node_id; 00218 double y1, y2, y3, sign_node_y1, sign_node_y2, sign_node_y3; 00219 double y1_term, y2_term, y3_term, y123_temp; 00220 00221 get_gauss_pts_and_weight(); 00222 00223 switch( numberNodes ) 00224 { 00225 case 8: 00226 for( i = 0; i < numberGaussPoints; i++ ) 00227 { 00228 for( j = 0; j < numberGaussPoints; j++ ) 00229 { 00230 for( k = 0; k < numberGaussPoints; k++ ) 00231 { 00232 y1 = gaussPointY[i]; 00233 y2 = gaussPointY[j]; 00234 y3 = gaussPointY[k]; 00235 00236 for( node_id = 0; node_id < numberNodes; node_id++ ) 00237 { 00238 get_signs_for_node_local_coord_hex( node_id, sign_node_y1, sign_node_y2, sign_node_y3 ); 00239 00240 y1_term = 1 + sign_node_y1 * y1; 00241 y2_term = 1 + sign_node_y2 * y2; 00242 y3_term = 1 + sign_node_y3 * y3; 00243 00244 shapeFunction[ife][node_id] = 0.125 * y1_term * y2_term * y3_term; 00245 dndy1GaussPts[ife][node_id] = 0.125 * sign_node_y1 * y2_term * y3_term; 00246 dndy2GaussPts[ife][node_id] = 0.125 * sign_node_y2 * y1_term * y3_term; 00247 dndy3GaussPts[ife][node_id] = 0.125 * sign_node_y3 * y1_term * y2_term; 00248 } 00249 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j] * gaussWeight[k]; 00250 ife++; 00251 } 00252 } 00253 } 00254 break; 00255 case 20: 00256 for( i = 0; i < numberGaussPoints; i++ ) 00257 { 00258 for( j = 0; j < numberGaussPoints; j++ ) 00259 { 00260 for( k = 0; k < numberGaussPoints; k++ ) 00261 { 00262 y1 = gaussPointY[i]; 00263 y2 = gaussPointY[j]; 00264 y3 = gaussPointY[k]; 00265 00266 for( node_id = 0; node_id < numberNodes; node_id++ ) 00267 { 00268 get_signs_for_node_local_coord_hex( node_id, sign_node_y1, sign_node_y2, sign_node_y3 ); 00269 00270 y1_term = 1 + sign_node_y1 * y1; 00271 y2_term = 1 + sign_node_y2 * y2; 00272 y3_term = 1 + sign_node_y3 * y3; 00273 y123_temp = sign_node_y1 * y1 + sign_node_y2 * y2 + sign_node_y3 * y3 - 2.; 00274 00275 switch( node_id ) 00276 { 00277 case 0: 00278 case 1: 00279 case 2: 00280 case 3: 00281 case 4: 00282 case 5: 00283 case 6: 00284 case 7: { 00285 shapeFunction[ife][node_id] = 0.125 * y1_term * y2_term * y3_term * y123_temp; 00286 dndy1GaussPts[ife][node_id] = 0.125 * sign_node_y1 * y123_temp * y2_term * y3_term + 00287 0.125 * y1_term * y2_term * y3_term * sign_node_y1; 00288 dndy2GaussPts[ife][node_id] = 0.125 * sign_node_y2 * y1_term * y3_term * y123_temp + 00289 0.125 * y1_term * y2_term * y3_term * sign_node_y2; 00290 dndy3GaussPts[ife][node_id] = 0.125 * sign_node_y3 * y1_term * y2_term * y123_temp + 00291 0.125 * y1_term * y2_term * y3_term * sign_node_y3; 00292 break; 00293 } 00294 case 8: 00295 case 10: 00296 case 16: 00297 case 18: { 00298 shapeFunction[ife][node_id] = 0.25 * ( 1 - y1 * y1 ) * y2_term * y3_term; 00299 dndy1GaussPts[ife][node_id] = -0.5 * y1 * y2_term * y3_term; 00300 dndy2GaussPts[ife][node_id] = 0.25 * ( 1 - y1 * y1 ) * sign_node_y2 * y3_term; 00301 dndy3GaussPts[ife][node_id] = 0.25 * ( 1 - y1 * y1 ) * y2_term * sign_node_y3; 00302 break; 00303 } 00304 case 9: 00305 case 11: 00306 case 17: 00307 case 19: { 00308 shapeFunction[ife][node_id] = 0.25 * ( 1 - y2 * y2 ) * y1_term * y3_term; 00309 dndy1GaussPts[ife][node_id] = 0.25 * ( 1 - y2 * y2 ) * sign_node_y1 * y3_term; 00310 dndy2GaussPts[ife][node_id] = -0.5 * y2 * y1_term * y3_term; 00311 dndy3GaussPts[ife][node_id] = 0.25 * ( 1 - y2 * y2 ) * y1_term * sign_node_y3; 00312 break; 00313 } 00314 case 12: 00315 case 13: 00316 case 14: 00317 case 15: { 00318 shapeFunction[ife][node_id] = 0.25 * ( 1 - y3 * y3 ) * y1_term * y2_term; 00319 dndy1GaussPts[ife][node_id] = 0.25 * ( 1 - y3 * y3 ) * sign_node_y1 * y2_term; 00320 dndy2GaussPts[ife][node_id] = 0.25 * ( 1 - y3 * y3 ) * y1_term * sign_node_y2; 00321 dndy3GaussPts[ife][node_id] = -0.5 * y3 * y1_term * y2_term; 00322 break; 00323 } 00324 } 00325 } 00326 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j] * gaussWeight[k]; 00327 ife++; 00328 } 00329 } 00330 } 00331 break; 00332 } 00333 } 00334 00335 void GaussIntegration::calculate_derivative_at_nodes( double dndy1_at_nodes[][maxNumberNodes], 00336 double dndy2_at_nodes[][maxNumberNodes] ) 00337 { 00338 double y1 = 0., y2 = 0.; 00339 int i; 00340 for( i = 0; i < numberNodes; i++ ) 00341 { 00342 switch( i ) 00343 { 00344 case 0: 00345 y1 = -1.; 00346 y2 = -1.; 00347 break; 00348 case 1: 00349 y1 = 1.; 00350 y2 = -1.; 00351 break; 00352 case 2: 00353 y1 = 1.; 00354 y2 = 1.; 00355 break; 00356 case 3: 00357 y1 = -1.; 00358 y2 = 1.; 00359 break; 00360 00361 // midside nodes if there is any 00362 00363 case 4: 00364 y1 = 0.; 00365 y2 = -1.; 00366 break; 00367 00368 case 5: 00369 y1 = 1.; 00370 y2 = 0.; 00371 break; 00372 00373 case 6: 00374 y1 = 0.; 00375 y2 = 1.; 00376 break; 00377 00378 case 7: 00379 y1 = -1.; 00380 y2 = 0.; 00381 break; 00382 } 00383 00384 switch( numberNodes ) 00385 { 00386 case 4: 00387 // dn_i/dy1 evaluated at node i 00388 dndy1_at_nodes[i][0] = -0.25 * ( 1 - y2 ); 00389 dndy1_at_nodes[i][1] = 0.25 * ( 1 - y2 ); 00390 dndy1_at_nodes[i][2] = 0.25 * ( 1 + y2 ); 00391 dndy1_at_nodes[i][3] = -0.25 * ( 1 + y2 ); 00392 00393 // dn_i/dy2 evaluated at node i 00394 dndy2_at_nodes[i][0] = -0.25 * ( 1 - y1 ); 00395 dndy2_at_nodes[i][1] = -0.25 * ( 1 + y1 ); 00396 dndy2_at_nodes[i][2] = 0.25 * ( 1 + y1 ); 00397 dndy2_at_nodes[i][3] = 0.25 * ( 1 - y1 ); 00398 break; 00399 00400 case 8: 00401 00402 dndy1_at_nodes[i][0] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 + y2 ); 00403 dndy1_at_nodes[i][1] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 - y2 ); 00404 dndy1_at_nodes[i][2] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 + y2 ); 00405 dndy1_at_nodes[i][3] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 - y2 ); 00406 00407 dndy1_at_nodes[i][4] = -y1 * ( 1 - y2 ); 00408 dndy1_at_nodes[i][5] = 0.5 * ( 1 - y2 * y2 ); 00409 dndy1_at_nodes[i][6] = -y1 * ( 1 + y2 ); 00410 dndy1_at_nodes[i][7] = -0.5 * ( 1 - y2 * y2 ); 00411 00412 dndy2_at_nodes[i][0] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 + y1 ); 00413 dndy2_at_nodes[i][1] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 - y1 ); 00414 dndy2_at_nodes[i][2] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 + y1 ); 00415 dndy2_at_nodes[i][3] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 - y1 ); 00416 00417 dndy2_at_nodes[i][4] = -0.5 * ( 1 - y1 * y1 ); 00418 dndy2_at_nodes[i][5] = -y2 * ( 1 + y1 ); 00419 dndy2_at_nodes[i][6] = 0.5 * ( 1 - y1 * y1 ); 00420 dndy2_at_nodes[i][7] = -y2 * ( 1 - y1 ); 00421 break; 00422 } 00423 } 00424 } 00425 00426 void GaussIntegration::calculate_derivative_at_nodes_3d( double dndy1_at_nodes[][maxNumberNodes], 00427 double dndy2_at_nodes[][maxNumberNodes], 00428 double dndy3_at_nodes[][maxNumberNodes] ) 00429 { 00430 double y1, y2, y3, sign_node_y1, sign_node_y2, sign_node_y3; 00431 double y1_term, y2_term, y3_term, y123_temp; 00432 int node_id, node_id_2; 00433 for( node_id = 0; node_id < numberNodes; node_id++ ) 00434 { 00435 get_signs_for_node_local_coord_hex( node_id, y1, y2, y3 ); 00436 00437 switch( numberNodes ) 00438 { 00439 case 8: 00440 for( node_id_2 = 0; node_id_2 < numberNodes; node_id_2++ ) 00441 { 00442 get_signs_for_node_local_coord_hex( node_id_2, sign_node_y1, sign_node_y2, sign_node_y3 ); 00443 y1_term = 1 + sign_node_y1 * y1; 00444 y2_term = 1 + sign_node_y2 * y2; 00445 y3_term = 1 + sign_node_y3 * y3; 00446 00447 dndy1_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y1 * y2_term * y3_term; 00448 00449 dndy2_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y2 * y1_term * y3_term; 00450 00451 dndy3_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y3 * y1_term * y2_term; 00452 } 00453 break; 00454 case 20: 00455 for( node_id_2 = 0; node_id_2 < numberNodes; node_id_2++ ) 00456 { 00457 get_signs_for_node_local_coord_hex( node_id_2, sign_node_y1, sign_node_y2, sign_node_y3 ); 00458 00459 y1_term = 1 + sign_node_y1 * y1; 00460 y2_term = 1 + sign_node_y2 * y2; 00461 y3_term = 1 + sign_node_y3 * y3; 00462 y123_temp = sign_node_y1 * y1 + sign_node_y2 * y2 + sign_node_y3 * y3 - 2.; 00463 switch( node_id_2 ) 00464 { 00465 case 0: 00466 case 1: 00467 case 2: 00468 case 3: 00469 case 4: 00470 case 5: 00471 case 6: 00472 case 7: { 00473 dndy1_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y1 * y2_term * y3_term * y123_temp + 00474 0.125 * y1_term * y2_term * y3_term * sign_node_y1; 00475 dndy2_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y2 * y1_term * y3_term * y123_temp + 00476 0.125 * y1_term * y2_term * y3_term * sign_node_y2; 00477 dndy3_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y3 * y1_term * y2_term * y123_temp + 00478 0.125 * y1_term * y2_term * y3_term * sign_node_y3; 00479 break; 00480 } 00481 case 8: 00482 case 10: 00483 case 16: 00484 case 18: { 00485 dndy1_at_nodes[node_id][node_id_2] = -0.5 * y1 * y2_term * y3_term; 00486 dndy2_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y1 * y1 ) * sign_node_y2 * y3_term; 00487 dndy3_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y1 * y1 ) * y2_term * sign_node_y3; 00488 break; 00489 } 00490 case 9: 00491 case 11: 00492 case 17: 00493 case 19: { 00494 dndy1_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y2 * y2 ) * sign_node_y1 * y3_term; 00495 dndy2_at_nodes[node_id][node_id_2] = -0.5 * y2 * y1_term * y3_term; 00496 dndy3_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y2 * y2 ) * y1_term * sign_node_y3; 00497 break; 00498 } 00499 case 12: 00500 case 13: 00501 case 14: 00502 case 15: { 00503 dndy1_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y3 * y3 ) * sign_node_y1 * y2_term; 00504 dndy2_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y3 * y3 ) * y1_term * sign_node_y2; 00505 dndy3_at_nodes[node_id][node_id_2] = -0.5 * y3 * y1_term * y2_term; 00506 break; 00507 } 00508 } 00509 } 00510 break; 00511 } 00512 } 00513 } 00514 00515 void GaussIntegration::get_signs_for_node_local_coord_hex( int node_id, double& sign_node_y1, double& sign_node_y2, 00516 double& sign_node_y3 ) 00517 { 00518 switch( node_id ) 00519 { 00520 case 0: 00521 sign_node_y1 = -1.; 00522 sign_node_y2 = -1.; 00523 sign_node_y3 = -1.; 00524 break; 00525 case 1: 00526 sign_node_y1 = 1.; 00527 sign_node_y2 = -1.; 00528 sign_node_y3 = -1.; 00529 break; 00530 case 2: 00531 sign_node_y1 = 1.; 00532 sign_node_y2 = 1.; 00533 sign_node_y3 = -1.; 00534 break; 00535 case 3: 00536 sign_node_y1 = -1.; 00537 sign_node_y2 = 1.; 00538 sign_node_y3 = -1.; 00539 break; 00540 case 4: 00541 sign_node_y1 = -1.; 00542 sign_node_y2 = -1.; 00543 sign_node_y3 = 1.; 00544 break; 00545 case 5: 00546 sign_node_y1 = 1.; 00547 sign_node_y2 = -1.; 00548 sign_node_y3 = 1.; 00549 break; 00550 case 6: 00551 sign_node_y1 = 1.; 00552 sign_node_y2 = 1.; 00553 sign_node_y3 = 1.; 00554 break; 00555 case 7: 00556 sign_node_y1 = -1.; 00557 sign_node_y2 = 1.; 00558 sign_node_y3 = 1.; 00559 break; 00560 case 8: 00561 sign_node_y1 = 0; 00562 sign_node_y2 = -1.; 00563 sign_node_y3 = -1.; 00564 break; 00565 case 9: 00566 sign_node_y1 = 1.; 00567 sign_node_y2 = 0; 00568 sign_node_y3 = -1.; 00569 break; 00570 case 10: 00571 sign_node_y1 = 0; 00572 sign_node_y2 = 1.; 00573 sign_node_y3 = -1.; 00574 break; 00575 case 11: 00576 sign_node_y1 = -1.; 00577 sign_node_y2 = 0.; 00578 sign_node_y3 = -1.; 00579 break; 00580 case 12: 00581 sign_node_y1 = -1.; 00582 sign_node_y2 = -1.; 00583 sign_node_y3 = 0.; 00584 break; 00585 case 13: 00586 sign_node_y1 = 1.; 00587 sign_node_y2 = -1.; 00588 sign_node_y3 = 0.; 00589 break; 00590 case 14: 00591 sign_node_y1 = 1.; 00592 sign_node_y2 = 1.; 00593 sign_node_y3 = 0.; 00594 break; 00595 case 15: 00596 sign_node_y1 = -1.; 00597 sign_node_y2 = 1.; 00598 sign_node_y3 = 0.; 00599 break; 00600 case 16: 00601 sign_node_y1 = 0; 00602 sign_node_y2 = -1.; 00603 sign_node_y3 = 1.; 00604 break; 00605 case 17: 00606 sign_node_y1 = 1.; 00607 sign_node_y2 = 0; 00608 sign_node_y3 = 1.; 00609 break; 00610 case 18: 00611 sign_node_y1 = 0; 00612 sign_node_y2 = 1.; 00613 sign_node_y3 = 1.; 00614 break; 00615 case 19: 00616 sign_node_y1 = -1.; 00617 sign_node_y2 = 0.; 00618 sign_node_y3 = 1.; 00619 break; 00620 } 00621 } 00622 00623 void GaussIntegration::get_tri_rule_pts_and_weight() 00624 { 00625 // get triangular rule integration points and weight 00626 00627 switch( numberGaussPoints ) 00628 { 00629 case 6: 00630 y1Area[0] = 0.09157621; 00631 y2Area[0] = 0.09157621; 00632 00633 y1Area[1] = 0.09157621; 00634 y2Area[1] = 0.8168476; 00635 00636 y1Area[2] = 0.8168476; 00637 y2Area[2] = 0.09157621; 00638 00639 y1Area[3] = 0.4459485; 00640 y2Area[3] = 0.4459485; 00641 00642 y1Area[4] = 0.4459485; 00643 y2Area[4] = 0.1081030; 00644 00645 y1Area[5] = 0.1081030; 00646 y2Area[5] = 0.4459485; 00647 00648 int i; 00649 for( i = 0; i < 3; i++ ) 00650 { 00651 totalGaussWeight[i] = 0.06348067; 00652 } 00653 for( i = 3; i < 6; i++ ) 00654 { 00655 totalGaussWeight[i] = 0.1289694; 00656 } 00657 break; 00658 } 00659 } 00660 00661 void GaussIntegration::calculate_shape_function_2d_tri() 00662 { 00663 int ife; 00664 double y1, y2, y3; 00665 get_tri_rule_pts_and_weight(); 00666 00667 for( ife = 0; ife < totalNumberGaussPts; ife++ ) 00668 { 00669 y1 = y1Area[ife]; 00670 y2 = y2Area[ife]; 00671 y3 = 1.0 - y1 - y2; 00672 00673 shapeFunction[ife][0] = y1 * ( 2. * y1 - 1. ); 00674 shapeFunction[ife][1] = y2 * ( 2. * y2 - 1. ); 00675 shapeFunction[ife][2] = y3 * ( 2. * y3 - 1. ); 00676 00677 shapeFunction[ife][3] = 4. * y1 * y2; 00678 shapeFunction[ife][4] = 4. * y2 * y3; 00679 shapeFunction[ife][5] = 4. * y1 * y3; 00680 00681 dndy1GaussPts[ife][0] = 4 * y1 - 1.; 00682 dndy1GaussPts[ife][1] = 0; 00683 dndy1GaussPts[ife][2] = 1 - 4. * y3; 00684 00685 dndy1GaussPts[ife][3] = 4. * y2; 00686 dndy1GaussPts[ife][4] = -4. * y2; 00687 dndy1GaussPts[ife][5] = 4. * ( 1 - 2 * y1 - y2 ); 00688 00689 dndy2GaussPts[ife][0] = 0.0; 00690 dndy2GaussPts[ife][1] = 4. * y2 - 1.; 00691 dndy2GaussPts[ife][2] = 1 - 4. * y3; 00692 00693 dndy2GaussPts[ife][3] = 4. * y1; 00694 dndy2GaussPts[ife][4] = 4. * ( 1 - y1 - 2. * y2 ); 00695 dndy2GaussPts[ife][5] = -4. * y1; 00696 } 00697 } 00698 00699 void GaussIntegration::calculate_derivative_at_nodes_2d_tri( double dndy1_at_nodes[][maxNumberNodes], 00700 double dndy2_at_nodes[][maxNumberNodes] ) 00701 { 00702 double y1 = 0., y2 = 0., y3; 00703 int i; 00704 for( i = 0; i < numberNodes; i++ ) 00705 { 00706 switch( i ) 00707 { 00708 case 0: 00709 y1 = 1.; 00710 y2 = 0.; 00711 break; 00712 case 1: 00713 y1 = 0.; 00714 y2 = 1.; 00715 break; 00716 case 2: 00717 y1 = 0.; 00718 y2 = 0.; 00719 break; 00720 case 3: 00721 y1 = 0.5; 00722 y2 = 0.5; 00723 break; 00724 case 4: 00725 y1 = 0.; 00726 y2 = 0.5; 00727 break; 00728 case 5: 00729 y1 = 0.5; 00730 y2 = 0.0; 00731 break; 00732 } 00733 00734 y3 = 1. - y1 - y2; 00735 00736 dndy1_at_nodes[i][0] = 4 * y1 - 1.; 00737 dndy1_at_nodes[i][1] = 0; 00738 dndy1_at_nodes[i][2] = 1 - 4. * y3; 00739 00740 dndy1_at_nodes[i][3] = 4. * y2; 00741 dndy1_at_nodes[i][4] = -4. * y2; 00742 dndy1_at_nodes[i][5] = 4. * ( 1 - 2 * y1 - y2 ); 00743 00744 dndy2_at_nodes[i][0] = 0.0; 00745 dndy2_at_nodes[i][1] = 4. * y2 - 1.; 00746 dndy2_at_nodes[i][2] = 1 - 4. * y3; 00747 00748 dndy2_at_nodes[i][3] = 4. * y1; 00749 dndy2_at_nodes[i][4] = 4. * ( 1 - y1 - 2. * y2 ); 00750 dndy2_at_nodes[i][5] = -4. * y1; 00751 } 00752 } 00753 void GaussIntegration::get_tet_rule_pts_and_weight() 00754 { 00755 // get tetrahedron rule integration points and weight 00756 00757 double a, b; 00758 switch( numberGaussPoints ) 00759 { 00760 case 1: 00761 // 1 integration point formula, degree of precision 1 00762 y1Volume[0] = 0.25; 00763 y2Volume[0] = 0.25; 00764 y3Volume[0] = 0.25; 00765 y4Volume[0] = 0.25; 00766 totalGaussWeight[0] = 1.; 00767 break; 00768 case 4: 00769 // 4 integration points formula, degree of precision 2 00770 a = 0.58541020; 00771 b = 0.13819660; 00772 00773 y1Volume[0] = a; 00774 y2Volume[0] = b; 00775 y3Volume[0] = b; 00776 y4Volume[0] = b; 00777 00778 y1Volume[1] = b; 00779 y2Volume[1] = a; 00780 y3Volume[1] = b; 00781 y4Volume[1] = b; 00782 00783 y1Volume[2] = b; 00784 y2Volume[2] = b; 00785 y3Volume[2] = a; 00786 y4Volume[2] = b; 00787 00788 y1Volume[3] = b; 00789 y2Volume[3] = b; 00790 y3Volume[3] = b; 00791 y4Volume[3] = a; 00792 00793 int i; 00794 for( i = 0; i < 4; i++ ) 00795 { 00796 totalGaussWeight[i] = 0.25; 00797 } 00798 break; 00799 } 00800 } 00801 00802 void GaussIntegration::calculate_shape_function_3d_tet() 00803 { 00804 int ife; 00805 double y1, y2, y3, y4; 00806 get_tet_rule_pts_and_weight(); 00807 00808 switch( numberNodes ) 00809 { 00810 case 10: // 10 nodes quadratic tet 00811 { 00812 for( ife = 0; ife < totalNumberGaussPts; ife++ ) 00813 { 00814 // y1,y2,y3,y4 are the volume coordinates 00815 y1 = y1Volume[ife]; 00816 y2 = y2Volume[ife]; 00817 y3 = y3Volume[ife]; 00818 y4 = y4Volume[ife]; 00819 00820 // shape function is the same as in ABAQUS 00821 // it is different from that in all the FEA book 00822 // in which node is the first node 00823 // here at node 1 y4=1 00824 shapeFunction[ife][0] = y4 * ( 2. * y4 - 1. ); 00825 shapeFunction[ife][1] = y1 * ( 2. * y1 - 1. ); 00826 shapeFunction[ife][2] = y2 * ( 2. * y2 - 1. ); 00827 shapeFunction[ife][3] = y3 * ( 2. * y3 - 1. ); 00828 00829 shapeFunction[ife][4] = 4. * y1 * y4; 00830 shapeFunction[ife][5] = 4. * y1 * y2; 00831 shapeFunction[ife][6] = 4. * y2 * y4; 00832 shapeFunction[ife][7] = 4. * y3 * y4; 00833 shapeFunction[ife][8] = 4. * y1 * y3; 00834 shapeFunction[ife][9] = 4. * y2 * y3; 00835 00836 dndy1GaussPts[ife][0] = 1 - 4 * y4; 00837 dndy1GaussPts[ife][1] = 4 * y1 - 1.; 00838 dndy1GaussPts[ife][2] = 0; 00839 dndy1GaussPts[ife][3] = 0; 00840 00841 dndy1GaussPts[ife][4] = 4. * ( y4 - y1 ); 00842 dndy1GaussPts[ife][5] = 4. * y2; 00843 dndy1GaussPts[ife][6] = -4. * y2; 00844 dndy1GaussPts[ife][7] = -4. * y3; 00845 dndy1GaussPts[ife][8] = 4. * y3; 00846 dndy1GaussPts[ife][9] = 0; 00847 00848 dndy2GaussPts[ife][0] = 1 - 4 * y4; 00849 dndy2GaussPts[ife][1] = 0; 00850 dndy2GaussPts[ife][2] = 4. * y2 - 1.; 00851 dndy2GaussPts[ife][3] = 0; 00852 00853 dndy2GaussPts[ife][4] = -4. * y1; 00854 dndy2GaussPts[ife][5] = 4. * y1; 00855 dndy2GaussPts[ife][6] = 4. * ( y4 - y2 ); 00856 dndy2GaussPts[ife][7] = -4. * y3; 00857 dndy2GaussPts[ife][8] = 0.; 00858 dndy2GaussPts[ife][9] = 4. * y3; 00859 00860 dndy3GaussPts[ife][0] = 1 - 4 * y4; 00861 dndy3GaussPts[ife][1] = 0; 00862 dndy3GaussPts[ife][2] = 0; 00863 dndy3GaussPts[ife][3] = 4. * y3 - 1.; 00864 00865 dndy3GaussPts[ife][4] = -4. * y1; 00866 dndy3GaussPts[ife][5] = 0; 00867 dndy3GaussPts[ife][6] = -4. * y2; 00868 dndy3GaussPts[ife][7] = 4. * ( y4 - y3 ); 00869 dndy3GaussPts[ife][8] = 4. * y1; 00870 dndy3GaussPts[ife][9] = 4. * y2; 00871 } 00872 break; 00873 } 00874 case 4: // four node linear tet for debug purpose 00875 { 00876 for( ife = 0; ife < totalNumberGaussPts; ife++ ) 00877 { 00878 y1 = y1Volume[ife]; 00879 y2 = y2Volume[ife]; 00880 y3 = y3Volume[ife]; 00881 y4 = y4Volume[ife]; 00882 00883 shapeFunction[ife][0] = y4; 00884 shapeFunction[ife][1] = y1; 00885 shapeFunction[ife][2] = y2; 00886 shapeFunction[ife][3] = y3; 00887 00888 dndy1GaussPts[ife][0] = -1.; 00889 dndy1GaussPts[ife][1] = 1; 00890 dndy1GaussPts[ife][2] = 0; 00891 dndy1GaussPts[ife][3] = 0; 00892 00893 dndy2GaussPts[ife][0] = -1.; 00894 dndy2GaussPts[ife][1] = 0; 00895 dndy2GaussPts[ife][2] = 1; 00896 dndy2GaussPts[ife][3] = 0; 00897 00898 dndy3GaussPts[ife][0] = -1.; 00899 dndy3GaussPts[ife][1] = 0; 00900 dndy3GaussPts[ife][2] = 0; 00901 dndy3GaussPts[ife][3] = 1; 00902 } 00903 break; 00904 } 00905 } 00906 } 00907 00908 void GaussIntegration::calculate_derivative_at_nodes_3d_tet( double dndy1_at_nodes[][maxNumberNodes], 00909 double dndy2_at_nodes[][maxNumberNodes], 00910 double dndy3_at_nodes[][maxNumberNodes] ) 00911 { 00912 double y1, y2, y3, y4; 00913 int i; 00914 00915 switch( numberNodes ) 00916 { 00917 case 10: { 00918 for( i = 0; i < numberNodes; i++ ) 00919 { 00920 get_node_local_coord_tet( i, y1, y2, y3, y4 ); 00921 00922 dndy1_at_nodes[i][0] = 1 - 4 * y4; 00923 dndy1_at_nodes[i][1] = 4 * y1 - 1.; 00924 dndy1_at_nodes[i][2] = 0; 00925 dndy1_at_nodes[i][3] = 0; 00926 00927 dndy1_at_nodes[i][4] = 4. * ( y4 - y1 ); 00928 dndy1_at_nodes[i][5] = 4. * y2; 00929 dndy1_at_nodes[i][6] = -4. * y2; 00930 dndy1_at_nodes[i][7] = -4. * y3; 00931 dndy1_at_nodes[i][8] = 4. * y3; 00932 dndy1_at_nodes[i][9] = 0; 00933 00934 dndy2_at_nodes[i][0] = 1 - 4 * y4; 00935 dndy2_at_nodes[i][1] = 0; 00936 dndy2_at_nodes[i][2] = 4. * y2 - 1.; 00937 dndy2_at_nodes[i][3] = 0; 00938 dndy2_at_nodes[i][4] = -4. * y1; 00939 dndy2_at_nodes[i][5] = 4. * y1; 00940 dndy2_at_nodes[i][6] = 4. * ( y4 - y2 ); 00941 dndy2_at_nodes[i][7] = -4. * y3; 00942 dndy2_at_nodes[i][8] = 0.; 00943 dndy2_at_nodes[i][9] = 4. * y3; 00944 00945 dndy3_at_nodes[i][0] = 1 - 4 * y4; 00946 dndy3_at_nodes[i][1] = 0; 00947 dndy3_at_nodes[i][2] = 0; 00948 dndy3_at_nodes[i][3] = 4. * y3 - 1.; 00949 00950 dndy3_at_nodes[i][4] = -4. * y1; 00951 dndy3_at_nodes[i][5] = 0; 00952 dndy3_at_nodes[i][6] = -4. * y2; 00953 dndy3_at_nodes[i][7] = 4. * ( y4 - y3 ); 00954 dndy3_at_nodes[i][8] = 4. * y1; 00955 dndy3_at_nodes[i][9] = 4. * y2; 00956 } 00957 break; 00958 } 00959 case 4: { 00960 for( i = 0; i < numberNodes; i++ ) 00961 { 00962 get_node_local_coord_tet( i, y1, y2, y3, y4 ); 00963 dndy1_at_nodes[i][0] = -1.; 00964 dndy1_at_nodes[i][1] = 1; 00965 dndy1_at_nodes[i][2] = 0; 00966 dndy1_at_nodes[i][3] = 0; 00967 00968 dndy2_at_nodes[i][0] = -1.; 00969 dndy2_at_nodes[i][1] = 0; 00970 dndy2_at_nodes[i][2] = 1; 00971 dndy2_at_nodes[i][3] = 0; 00972 00973 dndy3_at_nodes[i][0] = -1.; 00974 dndy3_at_nodes[i][1] = 0; 00975 dndy3_at_nodes[i][2] = 0; 00976 dndy3_at_nodes[i][3] = 1; 00977 } 00978 break; 00979 } 00980 } 00981 } 00982 00983 void GaussIntegration::get_node_local_coord_tet( int node_id, double& y1, double& y2, double& y3, double& y4 ) 00984 { 00985 switch( node_id ) 00986 { 00987 case 0: 00988 y1 = 0.; 00989 y2 = 0.; 00990 y3 = 0.; 00991 y4 = 1.; 00992 break; 00993 case 1: 00994 y1 = 1.; 00995 y2 = 0.; 00996 y3 = 0.; 00997 y4 = 0.; 00998 break; 00999 case 2: 01000 y1 = 0.; 01001 y2 = 1.; 01002 y3 = 0.; 01003 y4 = 0.; 01004 break; 01005 case 3: 01006 y1 = 0.; 01007 y2 = 0.; 01008 y3 = 1.; 01009 y4 = 0.; 01010 break; 01011 case 4: 01012 y1 = 0.5; 01013 y2 = 0.; 01014 y3 = 0.; 01015 y4 = 0.5; 01016 break; 01017 case 5: 01018 y1 = 0.5; 01019 y2 = 0.5; 01020 y3 = 0.; 01021 y4 = 0.; 01022 break; 01023 case 6: 01024 y1 = 0.; 01025 y2 = 0.5; 01026 y3 = 0.; 01027 y4 = 0.5; 01028 break; 01029 case 7: 01030 y1 = 0.; 01031 y2 = 0.0; 01032 y3 = 0.5; 01033 y4 = 0.5; 01034 break; 01035 case 8: 01036 y1 = 0.5; 01037 y2 = 0.; 01038 y3 = 0.5; 01039 y4 = 0.0; 01040 break; 01041 case 9: 01042 y1 = 0.; 01043 y2 = 0.5; 01044 y3 = 0.5; 01045 y4 = 0.; 01046 break; 01047 } 01048 }