MOAB: Mesh Oriented datABase  (version 5.2.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[], 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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines