![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 416
- Timestamp:
- 05/14/09 07:44:09 (16 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 2 edited
-
Icefront.cpp (modified) (5 diffs)
-
Icefront.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Icefront.cpp
r383 r416 255 255 if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ 256 256 CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type); 257 } 258 else if ((analysis_type==DiagnosticStokesAnalysisEnum())){ 259 260 CreatePVectorDiagnosticStokes( pg,inputs,analysis_type); 257 261 } 258 262 else{ … … 388 392 389 393 } 390 391 394 392 395 #undef __FUNCT__ … … 538 541 printf("pe_g->terms\n"); 539 542 for(i=0;i<numgridsload*NDOF2;i++){ 543 printf("%g ",*(pe_g->terms+i)); 544 } 545 } 546 #endif 547 548 /*Plug pe_g into vector: */ 549 VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); 550 551 /*Free ressources:*/ 552 xfree((void**)&element_nodes); 553 554 } 555 556 557 #undef __FUNCT__ 558 #define __FUNCT__ "Icefont::CreatePVectorDiagnosticStokes" 559 void Icefront::CreatePVectorDiagnosticStokes( Vec pg,void* vinputs, int analysis_type){ 560 561 int i,j; 562 563 const int numgrids=4; 564 const int NDOF4=4; 565 const int numdofs=numgrids*NDOF4; 566 int numberofdofspernode; 567 int doflist[numdofs]; 568 569 int fill; 570 double xyz_list[numgrids][3]; 571 double xyz_list_quad[numgrids+1][3]; //center node 572 573 double thickness_list[6]; //from penta element 574 double thickness_list_quad[6]; //selection of 4 thickness from penta elements 575 576 double bed_list[6]; //from penta element 577 double bed_list_quad[6]; //selection of 4 beds from penta elements 578 579 /*Objects: */ 580 double pe_g[numdofs]; 581 Matpar* matpar=NULL; 582 Node** element_nodes=NULL; 583 ParameterInputs* inputs=NULL; 584 585 /*quad grids: */ 586 int grid1,grid2,grid3,grid4; 587 588 /*normals: */ 589 double normal1[3]; 590 double normal2[3]; 591 double normal3[3]; 592 double normal4[3]; 593 double normal_norm; 594 double v15[3]; 595 double v25[3]; 596 double v35[3]; 597 double v45[3]; 598 599 inputs=(ParameterInputs*)vinputs; 600 601 /*check icefront is associated to a pentaelem: */ 602 if(element_type!=PentaEnum()){ 603 throw ErrorException(__FUNCT__," quad icefront is associated to a TriaElem!"); 604 } 605 /*Recover material and fill parameters: */ 606 matpar=(Matpar*)element->GetMatPar(); 607 fill=element->GetShelf(); 608 609 /* Set pe_g to 0: */ 610 for(i=0;i<numdofs;i++) pe_g[i]=0.0; 611 612 /* Get dof list and node coordinates: */ 613 GetDofList(&doflist[0],&numberofdofspernode); 614 GetElementNodeData( &xyz_list[0][0], nodes, numgrids); 615 616 //Identify which grids are comprised in the quad: 617 if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*)); 618 element->GetNodes((void**)element_nodes); 619 620 grid1=-1; grid2=-1; grid3=-1; grid4=-1; 621 for(i=0;i<6;i++){ 622 if (nodes[0]==element_nodes[i])grid1=i; 623 if (nodes[1]==element_nodes[i])grid2=i; 624 if (nodes[2]==element_nodes[i])grid3=i; 625 if (nodes[3]==element_nodes[i])grid4=i; 626 } 627 628 if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){ 629 throw ErrorException(__FUNCT__,"could not find element grids corresponding to quad icefront!"); 630 } 631 632 /*Build new xyz, bed and thickness lists for grids 1 to 4: */ 633 for(i=0;i<4;i++){ 634 for(j=0;j<3;j++){ 635 xyz_list_quad[i][j]=xyz_list[i][j]; 636 } 637 } 638 639 element->GetThicknessList(&thickness_list[0]); 640 element->GetBedList(&bed_list[0]); 641 642 thickness_list_quad[0]=thickness_list[grid1]; 643 thickness_list_quad[1]=thickness_list[grid2]; 644 thickness_list_quad[2]=thickness_list[grid3]; 645 thickness_list_quad[3]=thickness_list[grid4]; 646 647 bed_list_quad[0]=bed_list[grid1]; 648 bed_list_quad[1]=bed_list[grid2]; 649 bed_list_quad[2]=bed_list[grid3]; 650 bed_list_quad[3]=bed_list[grid4]; 651 652 //Create a new grid in the midle of the quad and add it at the end of the list 653 xyz_list_quad[4][0] = (xyz_list_quad[0][0]+xyz_list_quad[1][0]+xyz_list_quad[2][0]+xyz_list_quad[3][0])/4.0; 654 xyz_list_quad[4][1] = (xyz_list_quad[0][1]+xyz_list_quad[1][1]+xyz_list_quad[2][1]+xyz_list_quad[3][1])/4.0; 655 xyz_list_quad[4][2] = (xyz_list_quad[0][2]+xyz_list_quad[1][2]+xyz_list_quad[2][2]+xyz_list_quad[3][2])/4.0; 656 657 //Compute four normals (the quad is divided into four triangles) 658 v15[0]=xyz_list_quad[0][0]-xyz_list_quad[4][0]; 659 v15[1]=xyz_list_quad[0][1]-xyz_list_quad[4][1]; 660 v15[2]=xyz_list_quad[0][2]-xyz_list_quad[4][2]; 661 662 v25[0]=xyz_list_quad[1][0]-xyz_list_quad[4][0]; 663 v25[1]=xyz_list_quad[1][1]-xyz_list_quad[4][1]; 664 v25[2]=xyz_list_quad[1][2]-xyz_list_quad[4][2]; 665 666 v35[0]=xyz_list_quad[2][0]-xyz_list_quad[4][0]; 667 v35[1]=xyz_list_quad[2][1]-xyz_list_quad[4][1]; 668 v35[2]=xyz_list_quad[2][2]-xyz_list_quad[4][2]; 669 670 v45[0]=xyz_list_quad[3][0]-xyz_list_quad[4][0]; 671 v45[1]=xyz_list_quad[3][1]-xyz_list_quad[4][1]; 672 v45[2]=xyz_list_quad[3][2]-xyz_list_quad[4][2]; 673 674 cross(normal1,v15,v25); 675 cross(normal2,v25,v35); 676 cross(normal3,v35,v45); 677 cross(normal4,v45,v15); 678 679 normal_norm=norm(normal1);normal1[0]=normal1[0]/normal_norm;normal1[1]=normal1[1]/normal_norm;normal1[2]=normal1[2]/normal_norm; 680 normal_norm=norm(normal2);normal2[0]=normal2[0]/normal_norm;normal2[1]=normal2[1]/normal_norm;normal2[2]=normal2[2]/normal_norm; 681 normal_norm=norm(normal3);normal3[0]=normal3[0]/normal_norm;normal3[1]=normal3[1]/normal_norm;normal3[2]=normal3[2]/normal_norm; 682 normal_norm=norm(normal4);normal4[0]=normal4[0]/normal_norm;normal4[1]=normal4[1]/normal_norm;normal4[2]=normal4[2]/normal_norm; 683 684 685 //Compute load contribution for this quad: 686 QuadPressureLoadStokes(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0],fill); 687 688 689 #ifdef _DEBUGELEMENTS_ 690 if(my_rank==RANK && eid==ELID){ 691 printf("\nicefront load\n"); 692 printf("grids %i %i %i %i\n",grid1,grid2,grid3,grid4); 693 printf("rho_water %g\n",rho_water); 694 printf("rho_ice %g\n",rho_ice); 695 printf("gravity %g\n",gravity); 696 printf("thickness_list (%g,%g,%g,%g)\n",thickness_list[0],thickness_list[1],thickness_list[2],thickness_list[3]); 697 printf("bed_list (%g,%g,%g,%g)\n",bed_list[0],bed_list[1],bed_list[2],bed_list[3]); 698 printf("normal1 (%g,%g,%g)\n",normal1[0],normal1[1],normal1[2]); 699 printf("normal2 (%g,%g,%g)\n",normal2[0],normal2[1],normal2[2]); 700 printf("normal3 (%g,%g,%g)\n",normal3[0],normal3[1],normal3[2]); 701 printf("normal4 (%g,%g,%g)\n",normal4[0],normal4[1],normal4[2]); 702 printf("fill %i\n",fill); 703 printf("pe_g->terms\n"); 704 for(i=0;i<numgridsload*NDOF4;i++){ 540 705 printf("%g ",*(pe_g->terms+i)); 541 706 } … … 899 1064 complete_list[j][2]=0; 900 1065 } 901 tria->GetJacobianDeterminant 3d(&J[i],&complete_list[0][0],l1l2l3_tria);1066 tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],l1l2l3_tria); 902 1067 } 903 1068 … … 954 1119 delete tria; 955 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "Icefront::QuadPressureLoadStokes" 1124 void Icefront::QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 1125 double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill){ 1126 1127 1128 //The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows 1129 // 1130 // grid 4 +-----------------+ grid 3 1131 // |\2 1 /| 1132 // |1\ tria3 /2| 1133 // | \ / | 1134 // | \ / | 1135 // | \ / | 1136 // | \ / | 1137 // | \ 3 / | 1138 // |tria4 \ / 3 | 1139 // | 3 \grid5 | 1140 // | / \ | 1141 // | / 3 \ tria2| 1142 // | / \ | 1143 // | / \ | 1144 // | / \ | 1145 // | / tria1 \ | 1146 // |2/1 2\1| 1147 // grid1 +-----------------+ grid 2 1148 // 1149 // 1150 1151 int i,j; 1152 1153 double nx[4]; 1154 double ny[4]; 1155 double nz[4]; 1156 1157 /* gaussian points: */ 1158 int num_gauss,ig; 1159 double* first_gauss_area_coord = NULL; 1160 double* second_gauss_area_coord = NULL; 1161 double* third_gauss_area_coord = NULL; 1162 double* gauss_weights = NULL; 1163 double gauss_weight; 1164 double gauss_coord[3]; 1165 1166 double surface_list[5]; 1167 double x[5]; 1168 double y[5]; 1169 double z[5]; 1170 double l12,l14,l15,l23,l25,l34,l35,l45; 1171 double cos_theta_triangle1,cos_theta_triangle2,cos_theta_triangle3,cos_theta_triangle4; 1172 1173 double xyz_tria[12][2]; 1174 double l1l2l3_tria[3]; 1175 double r_tria,s_tria; 1176 double r_quad[4]; 1177 double s_quad[4]; 1178 double l1l4_tria[4][4]; 1179 1180 double J[4]; 1181 double z_g[4]; 1182 double ice_pressure_tria[4]; 1183 double air_pressure_tria; 1184 double water_level_above_g_tria; 1185 double water_pressure_tria; 1186 double pressure_tria[4]; 1187 1188 /*To use tria functionalities: */ 1189 Tria* tria=NULL; 1190 1191 /*Initialize fake tria: */ 1192 tria=new Tria(); 1193 1194 //Build the four normal vectors 1195 nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0]; 1196 ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1]; 1197 nz[0]=normal1[2]; nz[1]=normal2[2]; nz[2]=normal3[2]; nz[3]=normal4[2]; 1198 1199 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 1200 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 1201 1202 //Recover the surface of the four nodes 1203 for(i=0;i<4;i++){ 1204 surface_list[i]=thickness_list[i]+bed_list[i]; 1205 } 1206 //Add surface sor the fifth point (average of the surfaces) 1207 surface_list[4]=(surface_list[0]+surface_list[1]+surface_list[2]+surface_list[3])/4.0; 1208 1209 //Recover grid coordinates 1210 for(i=0;i<5;i++){ 1211 x[i]=*(xyz_list+3*i+0); 1212 y[i]=*(xyz_list+3*i+1); 1213 z[i]=*(xyz_list+3*i+2); 1214 } 1215 1216 //Build triangles in a 2D plan before using reference elements 1217 l12=pow((x[1]-x[0]),2)+pow((y[1]-y[0]),2)+pow((z[1]-z[0]),2); 1218 l14=pow((x[3]-x[0]),2)+pow((y[3]-y[0]),2)+pow((z[3]-z[0]),2); 1219 l15=pow((x[4]-x[0]),2)+pow((y[4]-y[0]),2)+pow((z[4]-z[0]),2); 1220 l23=pow((x[2]-x[1]),2)+pow((y[2]-y[1]),2)+pow((z[2]-z[1]),2); 1221 l25=pow((x[4]-x[1]),2)+pow((y[4]-y[1]),2)+pow((z[4]-z[1]),2); 1222 l34=pow((x[3]-x[2]),2)+pow((y[3]-y[2]),2)+pow((z[3]-z[2]),2); 1223 l35=pow((x[4]-x[2]),2)+pow((y[4]-y[2]),2)+pow((z[4]-z[2]),2); 1224 l45=pow((x[4]-x[3]),2)+pow((y[4]-y[3]),2)+pow((z[4]-z[3]),2); 1225 1226 cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15)); 1227 cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25)); 1228 cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35)); 1229 cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45)); 1230 1231 //First triangle : nodes 1, 2 and 5 1232 xyz_tria[0][0]=0; 1233 xyz_tria[0][1]=0; 1234 xyz_tria[1][0]=sqrt(l12); 1235 xyz_tria[1][1]=0; 1236 xyz_tria[2][0]=cos_theta_triangle1*sqrt(l15); 1237 xyz_tria[2][1]=sqrt(1-pow(cos_theta_triangle1,2))*sqrt(l15); 1238 1239 //Second triangle : nodes 2, 3 and 5 1240 xyz_tria[3][0]=0; 1241 xyz_tria[3][1]=0; 1242 xyz_tria[4][0]=sqrt(l23); 1243 xyz_tria[4][1]=0; 1244 xyz_tria[5][0]=cos_theta_triangle2*sqrt(l25); 1245 xyz_tria[5][1]=sqrt(1-pow(cos_theta_triangle2,2))*sqrt(l25); 1246 1247 //Third triangle : nodes 3, 4 and 5 1248 xyz_tria[6][0]=0; 1249 xyz_tria[6][1]=0; 1250 xyz_tria[7][0]=sqrt(l34); 1251 xyz_tria[7][1]=0; 1252 xyz_tria[8][0]=cos_theta_triangle3*sqrt(l35); 1253 xyz_tria[8][1]=sqrt(1-pow(cos_theta_triangle3,2))*sqrt(l35); 1254 1255 //Fourth triangle : nodes 4, 1 and 5 1256 xyz_tria[9][0]=0; 1257 xyz_tria[9][1]=0; 1258 xyz_tria[10][0]=sqrt(l14); 1259 xyz_tria[10][1]=0; 1260 xyz_tria[11][0]=cos_theta_triangle4*sqrt(l45); 1261 xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45); 1262 1263 1264 1265 //Start looping on the triangle's gaussian points: 1266 for(ig=0;ig<num_gauss;ig++){ 1267 1268 /*Pick up the gaussian point: */ 1269 gauss_weight=*(gauss_weights+ig); 1270 gauss_coord[0]=*(first_gauss_area_coord+ig); 1271 gauss_coord[1]=*(second_gauss_area_coord+ig); 1272 gauss_coord[2]=*(third_gauss_area_coord+ig); 1273 1274 tria->GetNodalFunctions(l1l2l3_tria,gauss_coord); 1275 1276 //Get the coordinates of gauss point for each triangle in penta/quad 1277 r_tria=gauss_coord[1]-gauss_coord[0]; 1278 s_tria=-3/sqrt(3)*(gauss_coord[0]+gauss_coord[1]-2/3); 1279 1280 //Coordinates of gauss points in the reference triangle 1281 r_quad[0]=r_tria; 1282 s_quad[0]=1/sqrt(3)*s_tria-2/3; 1283 r_quad[1]=-1/sqrt(3)*s_tria+2/3; 1284 s_quad[1]=r_tria; 1285 r_quad[2]=-r_tria; 1286 s_quad[2]=-1/sqrt(3)*s_tria+2/3; 1287 r_quad[3]=1/sqrt(3)*s_tria-2/3; 1288 s_quad[3]=-r_tria; 1289 1290 //Get the nodal function of the quad for the gauss points of each triangle 1291 for(i=0;i<4;i++){ 1292 l1l4_tria[i][0]=1.0/4.0*( 1293 (r_quad[i]-1.0)*(s_quad[i]-1.0) 1294 ); 1295 1296 l1l4_tria[i][1]=1.0/4.0*( 1297 -(r_quad[i]+1.0)*(s_quad[i]-1.0) 1298 ); 1299 1300 l1l4_tria[i][2]=1.0/4.0*( 1301 (r_quad[i]+1.0)*(s_quad[i]+1.0) 1302 ); 1303 1304 l1l4_tria[i][3]=1.0/4.0*( 1305 -(r_quad[i]-1.0)*(s_quad[i]+1.0) 1306 ); 1307 1308 } //for(i=0;i<4;i++) 1309 1310 1311 //Compute jacobian of each triangle 1312 for (i=0;i<4;i++){ 1313 double complete_list[3][3]; //a third coordinate is needed for the jacobian determinant calculation, here it is zero 1314 for(j=0;j<3;j++){ 1315 complete_list[j][0]=xyz_tria[3*i+j][0]; 1316 complete_list[j][1]=xyz_tria[3*i+j][1]; 1317 complete_list[j][2]=0; 1318 } 1319 tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],l1l2l3_tria); 1320 } 1321 1322 //Calculation of the z coordinate for the gaussian point ig for each triangle 1323 z_g[0]=z[0]*l1l2l3_tria[0]+z[1]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2]; 1324 z_g[1]=z[1]*l1l2l3_tria[0]+z[2]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2]; 1325 z_g[2]=z[2]*l1l2l3_tria[0]+z[3]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2]; 1326 z_g[3]=z[3]*l1l2l3_tria[0]+z[0]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2]; 1327 1328 //Loop on the triangles 1329 for(i=0;i<4;i++){ 1330 1331 //Loop on the grids of the quad 1332 //Calculate the ice pressure 1333 for(j=0;j<4;j++){ 1334 ice_pressure_tria[j]=gravity*rho_ice*(surface_list[j]-z_g[i]); 1335 } 1336 air_pressure_tria=0; 1337 1338 //Now deal with water pressure: 1339 if(fill==1){ //icefront ends in water 1340 water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level 1341 water_pressure_tria=rho_water*gravity*water_level_above_g_tria; 1342 } 1343 else if(fill==0){ 1344 water_pressure_tria=0; 1345 } 1346 else{ 1347 throw ErrorException(__FUNCT__,"QuadPressureLoadStokes error message: unknow fill type for icefront boundary condition"); 1348 } 1349 1350 //Add pressure from triangle i 1351 //Loop on the grids of the quad 1352 for(j=0;j<4;j++){ 1353 pressure_tria[j] = ice_pressure_tria[j] + water_pressure_tria + air_pressure_tria; 1354 } 1355 1356 1357 pe_g[0]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i]; 1358 pe_g[1]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i]; 1359 pe_g[2]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nz[i]; 1360 pe_g[3]+= 0; 1361 pe_g[4]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i]; 1362 pe_g[5]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i]; 1363 pe_g[6]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nz[i]; 1364 pe_g[7]+= 0; 1365 pe_g[8]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i]; 1366 pe_g[9]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i]; 1367 pe_g[10]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nz[i]; 1368 pe_g[11]+= 0; 1369 pe_g[12]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i]; 1370 pe_g[13]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i]; 1371 pe_g[14]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nz[i]; 1372 pe_g[15]+= 0; 1373 1374 1375 1376 } //for(i=0;i<4;i++) 1377 } //for(ig=0;ig<num_gauss;ig++) 1378 1379 /*Delete fake tria: */ 1380 delete tria; 1381 } -
issm/trunk/src/c/objects/Icefront.h
r246 r416 64 64 void CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type); 65 65 void CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type); 66 void CreatePVectorDiagnosticStokes( Vec pg,void* inputs, int analysis_type); 66 67 void GetDofList(int* doflist,int* pnumberofdofs); 67 68 void SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill); 68 69 void QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 70 double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill); 71 void QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 69 72 double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill); 70 73 void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
Note:
See TracChangeset
for help on using the changeset viewer.
