Changeset 5742
- Timestamp:
- 09/10/10 09:56:23 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Loads/Icefront.cpp
r5722 r5742 867 867 868 868 /* 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; 876 871 877 872 double surface_list[5]; … … 911 906 nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0]; 912 907 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);916 908 917 909 //Recover the surface of the four nodes … … 976 968 xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45); 977 969 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); 990 977 991 978 //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); 994 981 995 982 //Coordinates of gauss points in the reference triangle … … 1070 1057 1071 1058 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]; 1080 1067 1081 1068 1082 1069 1083 1070 } //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 1093 1073 delete tria; 1074 delete gauss; 1094 1075 } 1095 1076 /*}}}*/ … … 1129 1110 1130 1111 /* 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; 1138 1114 1139 1115 double surface_list[5]; … … 1172 1148 ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1]; 1173 1149 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);1177 1150 1178 1151 //Recover the surface of the four nodes … … 1237 1210 xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45); 1238 1211 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); 1251 1219 1252 1220 //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); 1255 1223 1256 1224 //Coordinates of gauss points in the reference triangle … … 1324 1292 pressure_tria = water_pressure_tria + air_pressure_tria; 1325 1293 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]; 1329 1297 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]; 1333 1301 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]; 1337 1305 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]; 1341 1309 pe_g[15]+= 0; 1342 1310 1343 1311 } //for(i=0;i<4;i++) 1344 1312 } //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 1353 1314 delete tria; 1315 delete gauss; 1354 1316 } 1355 1317 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.