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