Changeset 5720
- Timestamp:
- 09/09/10 12:58:29 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Loads/Icefront.cpp
r5719 r5720 504 504 void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){ 505 505 506 int i; 507 const int numnodes= NUMVERTICESSEG; 508 const int numdofs = numnodes *NDOF2; 509 int* doflist=NULL; 510 double xyz_list[numnodes][3]; 511 512 /*Objects: */ 513 double pe_g[numdofs] = {0.0}; 514 515 /*Segment*/ 516 double normal[2]; 517 int num_gauss; 518 double* segment_gauss_coord=NULL; 519 double* gauss_weights=NULL; 520 double gauss_weight; 521 double gauss_coord; 522 int ig; 523 double Jdet; 524 double L[2]; 525 double thickness,bed,pressure; 526 double ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity; 527 double surface_under_water,base_under_water; 528 int fill; 529 530 /*Inputs*/ 531 Input* thickness_input=NULL; 532 Input* bed_input=NULL; 533 Penta* penta=NULL; 534 535 BeamRef* beam=NULL; 536 537 //check that the element is onbed (collapsed formulation) otherwise:pe=0 538 if(!element->GetOnBed()) return; 539 540 /*Recover inputs: */ 541 thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum); 542 bed_input =((Penta*)element)->inputs->GetInput(BedEnum); 543 inputs->GetParameterValue(&fill,FillEnum); 544 545 /* Get dof list and node coordinates: */ 546 GetDofList(&doflist,MacAyealApproximationEnum); 547 GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes); 548 549 /*Get Matpar parameters*/ 550 rho_water=matpar->GetRhoWater(); 551 rho_ice =matpar->GetRhoIce(); 552 gravity =matpar->GetG(); 553 554 /*Get normal*/ 555 GetNormal(&normal[0],xyz_list); 556 557 //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions 558 num_gauss=3; 559 gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss); 560 561 for(ig=0;ig<num_gauss;ig++){ 562 563 /*Pick up the gaussian point: */ 564 gauss_weight=gauss_weights[ig]; 565 gauss_coord=segment_gauss_coord[ig]; 566 567 //compute Jacobian of segment 568 beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord); 569 570 /*Get thickness and bed*/ 571 this->GetParameterValue(&thickness,gauss_coord,thickness_input); 572 this->GetParameterValue(&bed, gauss_coord,bed_input); 573 574 /*Compute total pressure applied to ice front*/ 575 if (fill==WaterEnum){ 576 //icefront ends in water: 577 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 578 air_pressure=0; 579 580 //Now deal with water pressure 581 surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level 582 base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level 583 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); 584 } 585 else if (fill==AirEnum){ 586 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 587 air_pressure=0; 588 water_pressure=0; 589 } 590 else{ 591 ISSMERROR("fill type %i not supported yet",fill); 592 } 593 pressure = ice_pressure + water_pressure + air_pressure; 594 595 /*Get L*/ 596 beam->GetNodalFunctions(&L[0],gauss_coord); 597 598 for (i=0;i<numnodes;i++){ 599 pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i]; 600 pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i]; 601 } 602 } //for(ig=0;ig<num_gauss;ig++) 603 604 605 /*Plug pe_g into vector: */ 606 VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); 607 608 /*Free ressources:*/ 609 xfree((void**)&segment_gauss_coord); 610 xfree((void**)&gauss_weights); 611 xfree((void**)&doflist); 612 506 Icefront* icefront=NULL; 507 Penta* penta=NULL; 508 Tria* tria=NULL; 509 bool onbed; 510 511 /*Cast element onto Penta*/ 512 penta =(Penta*)this->element; 513 514 /*Return if not on bed*/ 515 penta->inputs->GetParameterValue(&onbed,ElementOnBedEnum); 516 if(!onbed) return; 517 518 /*Spawn Tria and call MacAyeal2d*/ 519 tria =(Tria*)penta->SpawnTria(0,1,2); 520 icefront=(Icefront*)this->copy(); 521 icefront->element=tria; 522 icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum)); 523 icefront->CreatePVectorDiagnosticMacAyeal2d(pg); 524 525 /*clean-up*/ 526 delete tria->matice; 527 delete tria; 528 delete icefront; 613 529 } 614 530 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.