Changeset 4946
- Timestamp:
- 08/03/10 10:31:37 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Loads
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Loads/Icefront.cpp
r4839 r4946 397 397 void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg){ 398 398 399 int i,j; 400 401 const int numgrids=2; 402 const int NDOF2=2; 403 int numberofdofspernode; 404 const int numdofs=numgrids*NDOF2; 405 int doflist[numdofs]; 406 double xyz_list[numgrids][3]; 407 408 double x1,y1,x2,y2; 409 410 /* input parameters: */ 411 double thickness_list_element[6]; 412 double thickness_list[2]; 413 double bed_list_element[6]; 414 double bed_list[2]; 415 416 int grid1,grid2; 417 double normal[2]; 418 double length; 419 int element_type; 399 int i; 400 const int numgrids = 2; 401 const int NDOF2 = 2; 402 const int numdofs = numgrids *NDOF2; 403 int numberofdofspernode; 404 int doflist[numdofs]; 405 double xyz_list[numgrids][3]; 420 406 421 407 /*Objects: */ 422 double pe_g[numdofs]={0.0}; 423 Matpar* matpar=NULL; 424 Node** element_nodes=NULL; 425 Node** nodes=NULL; 426 Element* element=NULL; 408 double pe_g[numdofs] = {0.0}; 409 Matpar *matpar = NULL; 410 Node **nodes = NULL; 411 Element *element = NULL; 412 413 /*Segment*/ 414 double normal[2]; 415 int num_gauss; 416 double* segment_gauss_coord=NULL; 417 double* gauss_weights=NULL; 418 double gauss_weight; 419 double gauss_coord; 420 int ig; 421 double Jdet; 422 double L[2]; 423 double thickness,bed,pressure; 424 double ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity; 425 double surface_under_water,base_under_water; 426 int fill; 427 428 /*Inputs*/ 429 Input* thickness_input=NULL; 430 Input* bed_input=NULL; 431 Tria* tria=NULL; 432 Penta* penta=NULL; 427 433 428 434 /*Recover hook objects: */ 429 matpar =(Matpar*)hmatpar->delivers();435 matpar =(Matpar*)hmatpar->delivers(); 430 436 element=(Element*)helement->delivers(); 431 nodes=(Node**)hnodes->deliverp(); 432 433 element_type=element->Enum(); 437 nodes =(Node**)hnodes->deliverp(); 438 BeamRef* beam=NULL; 434 439 435 440 //check that the element is onbed (collapsed formulation) otherwise:pe=0 436 if(element_type==PentaEnum){ 437 if (!element->GetOnBed()){ 438 return; 439 } 440 } 441 442 /*Identify which grids are comprised in the segment. */ 443 if(element_type==TriaEnum)element_nodes=(Node**)xmalloc(3*sizeof(Node*)); 444 if(element_type==PentaEnum)element_nodes=(Node**)xmalloc(6*sizeof(Node*)); 445 element->GetNodes((void**)element_nodes); 446 447 /*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */ 448 for(i=0;i<3;i++){ 449 if (nodes[0]==element_nodes[i])grid1=i; 450 if (nodes[1]==element_nodes[i])grid2=i; 451 } 452 441 if(element->Enum()==PentaEnum){ 442 if(!element->GetOnBed()) return; 443 thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum); 444 bed_input =((Penta*)element)->inputs->GetInput(BedEnum); 445 } 446 else{ 447 thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum); 448 bed_input =((Tria*)element)->inputs->GetInput(BedEnum); 449 } 450 451 /*Recover inputs: */ 452 inputs->GetParameterValue(&fill,FillEnum); 453 453 454 /* Get dof list and node coordinates: */ 454 455 GetDofList(&doflist[0],&numberofdofspernode); 455 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 456 457 /*Now build thickness_list_element and bed_list_element: */ 458 element->GetThicknessList(&thickness_list_element[0]); 459 element->GetBedList(&bed_list_element[0]); 460 461 /*Build thickness_list and bed_list: */ 462 thickness_list[0]=thickness_list_element[grid1]; 463 thickness_list[1]=thickness_list_element[grid2]; 464 bed_list[0]=bed_list_element[grid1]; 465 bed_list[1]=bed_list_element[grid2]; 466 467 //Recover grid coordinates 468 x1=xyz_list[0][0]; 469 y1=xyz_list[0][1]; 470 x2=xyz_list[1][0]; 471 y2=xyz_list[1][1]; 472 473 //Compute length and normal of segment 474 normal[0]=cos(atan2(x1-x2,y2-y1)); 475 normal[1]=sin(atan2(x1-x2,y2-y1)); 476 length=sqrt(pow(x2-x1,2)+pow(y2-y1,2)); 477 478 //Compute load contribution for this segment: 479 SegmentPressureLoad(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list,bed_list,normal,length); 456 GetVerticesCoordinates(&xyz_list[0][0],nodes,numgrids); 457 458 /*Get Matpar parameters*/ 459 rho_water=matpar->GetRhoWater(); 460 rho_ice=matpar->GetRhoIce(); 461 gravity=matpar->GetG(); 462 463 /*Get normal*/ 464 GetNormal(&normal[0],xyz_list); 465 466 //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions 467 num_gauss=3; 468 GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss); 469 470 for(ig=0;ig<num_gauss;ig++){ 471 472 /*Pick up the gaussian point: */ 473 gauss_weight=gauss_weights[ig]; 474 gauss_coord=segment_gauss_coord[ig]; 475 476 //compute Jacobian of segment 477 beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord); 478 479 /*Get thickness and bed*/ 480 this->GetParameterValue(&thickness,gauss_coord,thickness_input); 481 this->GetParameterValue(&bed, gauss_coord,bed_input); 482 483 /*Compute total pressure applied to ice front*/ 484 if (fill==WaterEnum){ 485 //icefront ends in water: 486 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 487 air_pressure=0; 488 489 //Now deal with water pressure 490 surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level 491 base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level 492 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); 493 } 494 else if (fill==AirEnum){ 495 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 496 air_pressure=0; 497 water_pressure=0; 498 } 499 else{ 500 ISSMERROR("fill type %i not supported yet",fill); 501 } 502 pressure = ice_pressure + water_pressure + air_pressure; 503 504 /*Get L*/ 505 beam->GetNodalFunctions(&L[0],gauss_coord); 506 507 for (i=0;i<numgrids;i++){ 508 pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i]; 509 pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i]; 510 } 511 } //for(ig=0;ig<num_gauss;ig++) 512 513 xfree((void**)&segment_gauss_coord); 514 xfree((void**)&gauss_weights); 480 515 481 516 /*Plug pe_g into vector: */ 482 517 VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); 483 484 /*Free ressources:*/485 xfree((void**)&element_nodes);486 518 487 519 } … … 1322 1354 } 1323 1355 /*}}}*/ 1324 /*FUNCTION Icefront::SegmentPressureLoad {{{1*/ 1325 1326 void Icefront::SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length){ 1327 1328 double nx,ny; 1329 double h1,h2,b1,b2; 1330 int num_gauss; 1331 double* segment_gauss_coord=NULL; 1332 double* gauss_weights=NULL; 1333 int ig; 1334 double Jdet; 1335 double thickness,bed; 1336 double ice_pressure,water_pressure,air_pressure; 1337 double surface_under_water,base_under_water; 1338 double pressure; 1339 int fill; 1340 1341 /*Recover inputs: */ 1342 inputs->GetParameterValue(&fill,FillEnum); 1343 1344 nx=normal[0]; 1345 ny=normal[1]; 1346 1347 1348 //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions 1349 num_gauss=3; 1350 GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss); 1351 1352 //recover thickness and bed at two grids 1353 h1=thickness_list[0]; 1354 h2=thickness_list[1]; 1355 b1=bed_list[0]; 1356 b2=bed_list[1]; 1357 1358 //compute Jacobian of segment 1359 Jdet=1./2.*length; 1360 1361 for(ig=0;ig<num_gauss;ig++){ 1362 1363 thickness=h1*(1+segment_gauss_coord[ig])/2+h2*(1-segment_gauss_coord[ig])/2; 1364 bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2; 1365 1366 if (fill==WaterEnum){ 1367 //icefront ends in water: 1368 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 1369 air_pressure=0; 1370 1371 //Now deal with water pressure 1372 surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level 1373 base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level 1374 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); 1375 } 1376 else if (fill==AirEnum){ 1377 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 1378 air_pressure=0; 1379 water_pressure=0; 1380 } 1381 else{ 1382 ISSMERROR("fill type %i not supported yet",fill); 1383 } 1384 1385 pressure = ice_pressure + water_pressure + air_pressure; 1386 1387 pe_g[2*0+0]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*nx; 1388 pe_g[2*0+1]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*ny; 1389 pe_g[2*1+0]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*nx; 1390 pe_g[2*1+1]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*ny; 1391 1392 } //for(ig=0;ig<num_gauss;ig++) 1393 1394 xfree((void**)&segment_gauss_coord); 1395 xfree((void**)&gauss_weights); 1396 } 1397 /*}}}*/ 1356 /*FUNCTION Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/ 1357 void Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){ 1358 1359 /*output: */ 1360 double value; 1361 1362 /*dynamic objects pointed to by hooks: */ 1363 Element* element=NULL; 1364 Node** nodes=NULL; 1365 1366 /*recover objects from hooks: */ 1367 element=(Element*)helement->delivers(); 1368 nodes=(Node**)hnodes->deliverp(); 1369 1370 /*Get value on Element 1*/ 1371 element->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in); 1372 1373 /*Assign output pointers:*/ 1374 *pvalue=value; 1375 } 1376 /*}}}*/ 1377 /*FUNCTION Icefront::GetNormal {{{1*/ 1378 void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){ 1379 1380 /*Build unit outward pointing vector*/ 1381 const int numgrids=2; 1382 double vector[2]; 1383 double norm; 1384 1385 vector[0]=xyz_list[1][0] - xyz_list[0][0]; 1386 vector[1]=xyz_list[1][1] - xyz_list[0][1]; 1387 1388 norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0)); 1389 1390 normal[0]= + vector[1]/norm; 1391 normal[1]= - vector[0]/norm; 1392 } 1393 /*}}}*/ -
issm/trunk/src/c/objects/Loads/Icefront.h
r4575 r4946 73 73 void CreatePVectorDiagnosticStokes( Vec pg); 74 74 void GetDofList(int* doflist,int* pnumberofdofs); 75 void SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length);76 75 void QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 77 76 double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list); 78 77 void QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 79 78 double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list); 79 void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in); 80 void GetNormal(double* normal,double xyz_list[2][3]); 80 81 /*}}}*/ 81 82 };
Note:
See TracChangeset
for help on using the changeset viewer.