Changeset 5742


Ignore:
Timestamp:
09/10/10 09:56:23 (15 years ago)
Author:
Mathieu Morlighem
Message:

removed old gauss

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5722 r5742  
    867867
    868868        /* gaussian points: */
    869         int     num_gauss,ig;
    870         double* first_gauss_area_coord  =  NULL;
    871         double* second_gauss_area_coord =  NULL;
    872         double* third_gauss_area_coord  =  NULL;
    873         double* gauss_weights           =  NULL;
    874         double  gauss_weight;
    875         double  gauss_coord[3];
     869        int     ig;
     870        GaussTria *gauss=NULL;
    876871
    877872        double  surface_list[5];
    … …  
    911906        nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0];
    912907        ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1];
    913 
    914         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    915         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    916908       
    917909        //Recover the surface of the four nodes
    … …  
    976968        xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45);
    977969
    978        
    979 
    980         //Start  looping on the triangle's gaussian points:
    981         for(ig=0;ig<num_gauss;ig++){
    982 
    983                 /*Pick up the gaussian point: */
    984                 gauss_weight=*(gauss_weights+ig);
    985                 gauss_coord[0]=*(first_gauss_area_coord+ig);
    986                 gauss_coord[1]=*(second_gauss_area_coord+ig);
    987                 gauss_coord[2]=*(third_gauss_area_coord+ig);
    988        
    989                 tria->GetNodalFunctions(l1l2l3_tria,gauss_coord);
     970        /* Start  looping on the number of gaussian points: */
     971        gauss=new GaussTria(2);
     972        for(ig=gauss->begin();ig<gauss->end();ig++){
     973
     974                gauss->GaussPoint(ig);
     975
     976                tria->GetNodalFunctions(l1l2l3_tria,gauss);
    990977
    991978                //Get the coordinates of gauss point for each triangle in penta/quad
    992                 r_tria=gauss_coord[1]-gauss_coord[0];
    993                 s_tria=-3/sqrt((double)3)*(gauss_coord[0]+gauss_coord[1]-2/3);
     979                r_tria=gauss->coord2-gauss->coord1;
     980                s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3);
    994981
    995982                //Coordinates of gauss points in the reference triangle
    … …  
    10701057
    10711058
    1072                         pe_g[0]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
    1073                         pe_g[1]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
    1074                         pe_g[2]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
    1075                         pe_g[3]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
    1076                         pe_g[4]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
    1077                         pe_g[5]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
    1078                         pe_g[6]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
    1079                         pe_g[7]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
     1059                        pe_g[0]+= J[i]*gauss->weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
     1060                        pe_g[1]+= J[i]*gauss->weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
     1061                        pe_g[2]+= J[i]*gauss->weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
     1062                        pe_g[3]+= J[i]*gauss->weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
     1063                        pe_g[4]+= J[i]*gauss->weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
     1064                        pe_g[5]+= J[i]*gauss->weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
     1065                        pe_g[6]+= J[i]*gauss->weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
     1066                        pe_g[7]+= J[i]*gauss->weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
    10801067
    10811068                       
    10821069
    10831070                } //for(i=0;i<4;i++)
    1084         } //for(ig=0;ig<num_gauss;ig++)
    1085        
    1086         /*Free ressources: */
    1087         xfree((void**)&first_gauss_area_coord);
    1088         xfree((void**)&second_gauss_area_coord);
    1089         xfree((void**)&third_gauss_area_coord);
    1090         xfree((void**)&gauss_weights);
    1091 
    1092         /*Delete fake tria: */
     1071        }
     1072       
    10931073        delete tria;
     1074        delete gauss;
    10941075}
    10951076/*}}}*/
    … …  
    11291110
    11301111        /* gaussian points: */
    1131         int     num_gauss,ig;
    1132         double* first_gauss_area_coord  =  NULL;
    1133         double* second_gauss_area_coord =  NULL;
    1134         double* third_gauss_area_coord  =  NULL;
    1135         double* gauss_weights           =  NULL;
    1136         double  gauss_weight;
    1137         double  gauss_coord[3];
     1112        int     ig;
     1113        GaussTria* gauss=NULL;
    11381114
    11391115        double  surface_list[5];
    … …  
    11721148        ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1];
    11731149        nz[0]=normal1[2]; nz[1]=normal2[2]; nz[2]=normal3[2]; nz[3]=normal4[2];
    1174 
    1175         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1176         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    11771150       
    11781151        //Recover the surface of the four nodes
    … …  
    12371210        xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45);
    12381211
    1239        
    1240 
    1241         //Start  looping on the triangle's gaussian points:
    1242         for(ig=0;ig<num_gauss;ig++){
    1243 
    1244                 /*Pick up the gaussian point: */
    1245                 gauss_weight=*(gauss_weights+ig);
    1246                 gauss_coord[0]=*(first_gauss_area_coord+ig);
    1247                 gauss_coord[1]=*(second_gauss_area_coord+ig);
    1248                 gauss_coord[2]=*(third_gauss_area_coord+ig);
    1249        
    1250                 tria->GetNodalFunctions(l1l2l3_tria,gauss_coord);
     1212        /* Start  looping on the number of gaussian points: */
     1213        gauss=new GaussTria(2);
     1214        for(ig=gauss->begin();ig<gauss->end();ig++){
     1215
     1216                gauss->GaussPoint(ig);
     1217
     1218                tria->GetNodalFunctions(l1l2l3_tria,gauss);
    12511219
    12521220                //Get the coordinates of gauss point for each triangle in penta/quad
    1253                 r_tria=gauss_coord[1]-gauss_coord[0];
    1254                 s_tria=-3/sqrt((double)3)*(gauss_coord[0]+gauss_coord[1]-2/3);
     1221                r_tria=gauss->coord2-gauss->coord1;
     1222                s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3);
    12551223
    12561224                //Coordinates of gauss points in the reference triangle
    … …  
    13241292                        pressure_tria = water_pressure_tria + air_pressure_tria;
    13251293
    1326                         pe_g[0]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nx[i];
    1327                         pe_g[1]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*ny[i];
    1328                         pe_g[2]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nz[i];
     1294                        pe_g[0]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nx[i];
     1295                        pe_g[1]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*ny[i];
     1296                        pe_g[2]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nz[i];
    13291297                        pe_g[3]+= 0;
    1330                         pe_g[4]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nx[i];
    1331                         pe_g[5]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*ny[i];
    1332                         pe_g[6]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nz[i];
     1298                        pe_g[4]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nx[i];
     1299                        pe_g[5]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*ny[i];
     1300                        pe_g[6]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nz[i];
    13331301                        pe_g[7]+= 0;
    1334                         pe_g[8]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nx[i];
    1335                         pe_g[9]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*ny[i];
    1336                         pe_g[10]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nz[i];
     1302                        pe_g[8]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nx[i];
     1303                        pe_g[9]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*ny[i];
     1304                        pe_g[10]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nz[i];
    13371305                        pe_g[11]+= 0;
    1338                         pe_g[12]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nx[i];
    1339                         pe_g[13]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*ny[i];
    1340                         pe_g[14]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nz[i];
     1306                        pe_g[12]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nx[i];
     1307                        pe_g[13]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*ny[i];
     1308                        pe_g[14]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nz[i];
    13411309                        pe_g[15]+= 0;
    13421310
    13431311                } //for(i=0;i<4;i++)
    13441312        } //for(ig=0;ig<num_gauss;ig++)
    1345        
    1346         /*Free ressources: */
    1347         xfree((void**)&first_gauss_area_coord);
    1348         xfree((void**)&second_gauss_area_coord);
    1349         xfree((void**)&third_gauss_area_coord);
    1350         xfree((void**)&gauss_weights);
    1351 
    1352         /*Delete fake tria: */
     1313
    13531314        delete tria;
     1315        delete gauss;
    13541316}
    13551317/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.