Changeset 5715
- Timestamp:
- 09/09/10 09:23:48 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
r5707 r5715 115 115 RiftfrontEnum, 116 116 SegmentRiftfrontEnum, 117 MacAyealIceFrontEnum, 117 MacAyeal2dIceFrontEnum, 118 MacAyeal3dIceFrontEnum, 118 119 PattynIceFrontEnum, 119 120 StokesIceFrontEnum, -
issm/trunk/src/c/EnumDefinitions/EnumToString.cpp
r5707 r5715 100 100 case RiftfrontEnum : return "Riftfront"; 101 101 case SegmentRiftfrontEnum : return "SegmentRiftfront"; 102 case MacAyealIceFrontEnum : return "MacAyealIceFront"; 102 case MacAyeal2dIceFrontEnum : return "MacAyeal2dIceFront"; 103 case MacAyeal3dIceFrontEnum : return "MacAyeal3dIceFront"; 103 104 case PattynIceFrontEnum : return "PattynIceFront"; 104 105 case StokesIceFrontEnum : return "StokesIceFront"; -
issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp
r5707 r5715 98 98 else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; 99 99 else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; 100 else if (strcmp(name,"MacAyealIceFront")==0) return MacAyealIceFrontEnum; 100 else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum; 101 else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum; 101 102 else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum; 102 103 else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum; -
issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
r5607 r5715 60 60 61 61 /*Create and add load: */ 62 if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum)){ 63 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyealIceFrontEnum,DiagnosticHorizAnalysisEnum)); 62 if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==2){ 63 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum)); 64 count++; 65 } 66 else if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==3){ 67 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum)); 64 68 count++; 65 69 } … … 73 77 } 74 78 else if ((int)*(iomodel->elements_type+element)==(MacAyealPattynApproximationEnum)){ 75 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal IceFrontEnum,DiagnosticHorizAnalysisEnum));79 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum)); 76 80 count++; 77 81 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum)); -
issm/trunk/src/c/objects/Loads/Icefront.cpp
r5714 r5715 68 68 icefront_mparid=iomodel->numberofelements+1; //matlab indexing 69 69 70 if (in_icefront_type==MacAyeal IceFrontEnum){70 if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){ 71 71 icefront_node_ids[0]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+0); 72 72 icefront_node_ids[1]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+1); … … 330 330 int type; 331 331 inputs->GetParameterValue(&type,TypeEnum); 332 if (type==MacAyealIceFrontEnum){ 333 CreatePVectorDiagnosticMacAyeal( pg); 332 if (type==MacAyeal2dIceFrontEnum){ 333 CreatePVectorDiagnosticMacAyeal2d( pg); 334 } 335 else if (type==MacAyeal3dIceFrontEnum){ 336 CreatePVectorDiagnosticMacAyeal3d( pg); 334 337 } 335 338 else if (type==PattynIceFrontEnum){ … … 419 422 420 423 /*Icefront numerics: */ 421 /*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal {{{1*/422 void Icefront::CreatePVectorDiagnosticMacAyeal ( Vec pg){424 /*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/ 425 void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){ 423 426 424 427 int i; … … 450 453 Input* bed_input=NULL; 451 454 Tria* tria=NULL; 452 Penta* penta=NULL;453 455 454 456 BeamRef* beam=NULL; 455 457 456 //check that the element is onbed (collapsed formulation) otherwise:pe=0457 if(element->Enum()==PentaEnum){458 if(!element->GetOnBed()) return;459 thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);460 bed_input =((Penta*)element)->inputs->GetInput(BedEnum);461 }462 else{463 thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);464 bed_input =((Tria*)element)->inputs->GetInput(BedEnum);465 }466 467 458 /*Recover inputs: */ 459 thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum); 460 bed_input =((Tria*)element)->inputs->GetInput(BedEnum); 468 461 inputs->GetParameterValue(&fill,FillEnum); 469 462 … … 514 507 } 515 508 else{ 516 ISSMERROR("fill type % i not supported yet",fill);509 ISSMERROR("fill type %s not supported yet",EnumToString(fill)); 517 510 } 518 511 pressure = ice_pressure + water_pressure + air_pressure; … … 528 521 529 522 523 /*Plug pe_g into vector: */ 524 VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); 525 526 /*Free ressources:*/ 527 xfree((void**)&segment_gauss_coord); 528 xfree((void**)&gauss_weights); 529 xfree((void**)&doflist); 530 531 } 532 /*}}}*/ 533 /*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{1*/ 534 void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){ 535 536 int i; 537 const int numnodes= NUMVERTICESSEG; 538 const int numdofs = numnodes *NDOF2; 539 int* doflist=NULL; 540 double xyz_list[numnodes][3]; 541 542 /*Objects: */ 543 double pe_g[numdofs] = {0.0}; 544 545 /*Segment*/ 546 double normal[2]; 547 int num_gauss; 548 double* segment_gauss_coord=NULL; 549 double* gauss_weights=NULL; 550 double gauss_weight; 551 double gauss_coord; 552 int ig; 553 double Jdet; 554 double L[2]; 555 double thickness,bed,pressure; 556 double ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity; 557 double surface_under_water,base_under_water; 558 int fill; 559 560 /*Inputs*/ 561 Input* thickness_input=NULL; 562 Input* bed_input=NULL; 563 Penta* penta=NULL; 564 565 BeamRef* beam=NULL; 566 567 //check that the element is onbed (collapsed formulation) otherwise:pe=0 568 if(!element->GetOnBed()) return; 569 570 /*Recover inputs: */ 571 thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum); 572 bed_input =((Penta*)element)->inputs->GetInput(BedEnum); 573 inputs->GetParameterValue(&fill,FillEnum); 574 575 /* Get dof list and node coordinates: */ 576 GetDofList(&doflist,MacAyealApproximationEnum); 577 GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes); 578 579 /*Get Matpar parameters*/ 580 rho_water=matpar->GetRhoWater(); 581 rho_ice =matpar->GetRhoIce(); 582 gravity =matpar->GetG(); 583 584 /*Get normal*/ 585 GetNormal(&normal[0],xyz_list); 586 587 //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions 588 num_gauss=3; 589 gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss); 590 591 for(ig=0;ig<num_gauss;ig++){ 592 593 /*Pick up the gaussian point: */ 594 gauss_weight=gauss_weights[ig]; 595 gauss_coord=segment_gauss_coord[ig]; 596 597 //compute Jacobian of segment 598 beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord); 599 600 /*Get thickness and bed*/ 601 this->GetParameterValue(&thickness,gauss_coord,thickness_input); 602 this->GetParameterValue(&bed, gauss_coord,bed_input); 603 604 /*Compute total pressure applied to ice front*/ 605 if (fill==WaterEnum){ 606 //icefront ends in water: 607 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 608 air_pressure=0; 609 610 //Now deal with water pressure 611 surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level 612 base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level 613 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); 614 } 615 else if (fill==AirEnum){ 616 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 617 air_pressure=0; 618 water_pressure=0; 619 } 620 else{ 621 ISSMERROR("fill type %i not supported yet",fill); 622 } 623 pressure = ice_pressure + water_pressure + air_pressure; 624 625 /*Get L*/ 626 beam->GetNodalFunctions(&L[0],gauss_coord); 627 628 for (i=0;i<numnodes;i++){ 629 pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i]; 630 pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i]; 631 } 632 } //for(ig=0;ig<num_gauss;ig++) 633 634 530 635 /*Plug pe_g into vector: */ 531 636 VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); … … 818 923 819 924 /*How many nodes? :*/ 820 if(type==MacAyealIceFrontEnum)numberofnodes=2; 821 else numberofnodes=4; 925 if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum) 926 numberofnodes=2; 927 else 928 numberofnodes=4; 822 929 823 930 /*Figure out size of doflist: */ -
issm/trunk/src/c/objects/Loads/Icefront.h
r5714 r5715 76 76 /*}}}*/ 77 77 /*Load management: {{{1*/ 78 void CreatePVectorDiagnosticMacAyeal( Vec pg); 79 void CreatePVectorDiagnosticPattyn( Vec pg); 80 void CreatePVectorDiagnosticStokes( Vec pg); 78 void CreatePVectorDiagnosticMacAyeal2d(Vec pg); 79 void CreatePVectorDiagnosticMacAyeal3d(Vec pg); 80 void CreatePVectorDiagnosticPattyn(Vec pg); 81 void CreatePVectorDiagnosticStokes(Vec pg); 81 82 void GetDofList(int** pdoflist,int approximation_enum=0); 82 83 void QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
Note:
See TracChangeset
for help on using the changeset viewer.