Index: /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp (revision 7322) +++ /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp (revision 7323) @@ -24,5 +24,5 @@ /*call different migration modules, according to user wishes: */ if(migration_style==AgressiveMigrationEnum) AgressiveMigration(elements,nodes, vertices,loads,materials, parameters); - else if(migration_style==SoftMigrationEnum) AgressiveMigration(elements,nodes, vertices,loads,materials, parameters); + else if(migration_style==SoftMigrationEnum) SoftMigration(elements,nodes, vertices,loads,materials, parameters); else if(migration_style==NoneEnum)_printf_("%s\n","NoneEnum supplied for migration style, doing nothing!"); else _error_("%s not supported yet!",EnumToString(migration_style)); Index: /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h =================================================================== --- /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h (revision 7322) +++ /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h (revision 7323) @@ -11,12 +11,12 @@ /* local prototypes: */ - -void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); -void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); - -double* CreateElementOnIceShelf(Elements* elements); -bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf); -double* CreateElementTouchingIceShelf(Elements* elements); -int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf); - +void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); +void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); +double* PotentialSheetUngrounding(Elements* elements,Nodes* nodes,Parameters* parameters); +double* PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding); +bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf); +double* CreateElementOnIceShelf(Elements* elements); +double* CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf); +int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf); +Vec CreateNodesOnIceShelf(Nodes* nodes,int analysis_type); #endif /* _GROUNDINGLINEMIGRATIONXLOCAL_H */ Index: /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp (revision 7322) +++ /issm/trunk/src/c/modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp (revision 7323) @@ -10,5 +10,5 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ +void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ //{{{1 /*Here, whatever grids inside the ice sheet want to unground, we allow -> instantaneous transmission of water through the bedrock: */ @@ -28,10 +28,141 @@ for(i=0;iSize();i++){ element=(Element*)elements->GetObjectByOffset(i); - element->AgressiveShelfSync(); - } - -} - -bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf){ + element->ShelfSync(); + } + +} +/*}}}*/ + +void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ //{{{1 + + + /*intermediary: */ + int i; + double* potential_sheet_ungrounding=NULL; + double* sheet_ungrounding=NULL; + Element* element=NULL; + + + _printf_(VerboseModule()," Migrating grounding line\n"); + + /*First, figure out which nodes that are on the ice sheet want to unground. Build a flags vector for this (size number of nodes: */ + potential_sheet_ungrounding=PotentialSheetUngrounding(elements,nodes,parameters); + + /*Now, propagate ice shelf into connex areas of the ice sheet that potentially want to unground: */ + sheet_ungrounding=PropagateShelfIntoConnexIceSheet(elements,nodes,parameters,potential_sheet_ungrounding); + + /*Now, use the sheet_ungrounding flags to unground the ice sheet (at the same time, take care of grounding elements of the ice shelf + * that want to: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->SoftMigration(sheet_ungrounding); + } + + /*Synchronize shelf status: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->ShelfSync(); + } + + /*free ressouces: */ + xfree((void**)&potential_sheet_ungrounding); + xfree((void**)&sheet_ungrounding); + +} +/*}}}*/ +double* PotentialSheetUngrounding(Elements* elements,Nodes* nodes,Parameters* parameters){ //{{{1 + + int i; + Element* element=NULL; + Vec vec_potential_sheet_ungrounding=NULL; + double* potential_sheet_ungrounding=NULL; + int numnods; + int analysis_type; + + /*recover parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */ + numnods=nodes->NumberOfNodes(analysis_type); + vec_potential_sheet_ungrounding=NewVec(numnods); + + /*Loop through elements, and fill vec_potential_sheet_ungrounding: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->PotentialSheetUngrounding(vec_potential_sheet_ungrounding); + } + + /*Assemble vector: */ + VecAssemblyBegin(vec_potential_sheet_ungrounding); + VecAssemblyEnd(vec_potential_sheet_ungrounding); + + /*Serialize vector: */ + VecToMPISerial(&potential_sheet_ungrounding,vec_potential_sheet_ungrounding); + + /*free ressouces: */ + VecFree(&vec_potential_sheet_ungrounding); + + return potential_sheet_ungrounding; +} +/*}}}*/ +double* PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding){ //{{{1 + + int i; + Element* element=NULL; + int numnods; + int analysis_type; + Vec vec_nodes_on_iceshelf=NULL; + double* nodes_on_iceshelf=NULL; + int nflipped,local_nflipped; + double* elements_touching_iceshelf=NULL; + + /*recover parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*recover vec_nodes_on_iceshelf: */ + vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,analysis_type); + + nflipped=1; //bootstrap + while(nflipped){ + + /*get a list of potential elements that have nodes on ice shelf: */ + elements_touching_iceshelf=CreateElementTouchingIceShelf(elements,vec_nodes_on_iceshelf); + + /*now, go through elements_touching_iceshelf, and if they have nodes inside potential_sheet_ungrounding, + * flag it: */ + local_nflipped=0; + + /*serialize vec_nodes_on_iceshelf, needed by element's UpdatePotentialSheetUngrounding routine, to figure out if + * nodes have flipped from grounded to ungrounded: */ + VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf); + + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + if(elements_touching_iceshelf[element->Sid()]){ + nflipped+=element->UpdatePotentialSheetUngrounding(potential_sheet_ungrounding,vec_nodes_on_iceshelf,nodes_on_iceshelf); + } + } + + MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); + //_printf_(VerboseConvergence()," number of grounded vertices connected to grounding line: %i\n",nflipped); + _printf_(1," number of grounded vertices connected to grounding line: %i\n",nflipped); + + /*Avoid leaks: */ + xfree((void**)&nodes_on_iceshelf); + } + + /*Serialize: */ + VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf); + + /*Free ressources:*/ + VecFree(&vec_nodes_on_iceshelf); + xfree((void**)&elements_touching_iceshelf); + + return nodes_on_iceshelf; +} +/*}}}*/ + + +bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf){ //{{{1 int i; @@ -67,6 +198,6 @@ return element_on_gl; } - -double* CreateElementOnIceShelf(Elements* elements){ +/*}}}*/ +double* CreateElementOnIceShelf(Elements* elements){ //{{{1 int i; @@ -97,7 +228,6 @@ return element_on_iceshelf; } - - -double* CreateElementTouchingIceShelf(Elements* elements){ +/*}}}*/ +double* CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf){ //{{{1 int i; @@ -105,4 +235,8 @@ Vec vec_element_touching_iceshelf=NULL; double* element_touching_iceshelf=NULL; + double* nodes_on_iceshelf=NULL; + + /*serialize: */ + VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf); /*Create vector holding all the elements IsOnShelf flags: */ @@ -112,5 +246,5 @@ for(i=0;iSize();i++){ element=(Element*)elements->GetObjectByOffset(i); - VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelf()?1.0:0.0,INSERT_VALUES); + VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_iceshelf)?1.0:0.0,INSERT_VALUES); } @@ -124,11 +258,11 @@ /*free ressouces: */ VecFree(&vec_element_touching_iceshelf); + xfree((void**)&nodes_on_iceshelf); return element_touching_iceshelf; } - - -int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){ - +/*}}}*/ +int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){//{{{1 + int i; Element* element=NULL; @@ -175,49 +309,31 @@ return nflipped; } - - -void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ - - int i,j; - Element* element=NULL; - double* element_touching_iceshelf=NULL; - int nflipped; - - _printf_(VerboseModule()," Migrating grounding line\n"); - - /*Loop until no more nodes and elements change from grounded to floating: */ - nflipped=1; //arbitrary, get things started - - int count=0; - while(nflipped){ - - if (count==4)break; - - /*reset counter: */ - nflipped=0; - - /*Create vector holding all the elements that touch the ice shelf, by any node: */ - element_touching_iceshelf=CreateElementTouchingIceShelf(elements); - - - /*Carry out grounding line migration for those elements: */ - for(i=0;iSize();i++){ - element=(Element*)elements->GetObjectByOffset(i); - if(element_touching_iceshelf[element->Sid()]) element->MigrateGroundingLine(); - } - - /*Now, update shelf flags in nodes and elements: */ - nflipped=UpdateShelfStatus(elements,nodes,parameters,element_touching_iceshelf); - //_printf_(VerboseModule()," number of migrated nodes: %i\n",nflipped); - extern int my_rank;if(my_rank==0)printf(" number of migrated nodes: %i\n",nflipped); - - /*avoid memory leaks: */ - xfree((void**)&element_touching_iceshelf); - count++; - } - - - /*free ressouces: */ - xfree((void**)&element_touching_iceshelf); - -} +/*}}}*/ +Vec CreateNodesOnIceShelf(Nodes* nodes,int analysis_type){ //{{{1 + + /*output: */ + Vec nodes_on_iceshelf=NULL; + + /*intermediary: */ + int numnods; + int i; + Node* node=NULL; + + + /*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */ + numnods=nodes->NumberOfNodes(analysis_type); + nodes_on_iceshelf=NewVec(numnods); + + /*Loop through nodes, and fill nodes_on_iceshelf: */ + for(i=0;iSize();i++){ + node=(Node*)nodes->GetObjectByOffset(i); + if(node->IsOnShelf())VecSetValue(nodes_on_iceshelf,node->Sid()+1,1,INSERT_VALUES); + } + + /*Assemble vector: */ + VecAssemblyBegin(nodes_on_iceshelf); + VecAssemblyEnd(nodes_on_iceshelf); + + return nodes_on_iceshelf; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 7322) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 7323) @@ -35,4 +35,5 @@ virtual bool IsOnShelf()=0; virtual bool IsNodeOnShelf()=0; + virtual bool IsNodeOnShelfFromFlags(double* flags)=0; virtual bool IsOnBed()=0; virtual void GetParameterListOnVertices(double* pvalue,int enumtype)=0; @@ -84,8 +85,11 @@ virtual double TimeAdapt()=0; virtual void AgressiveMigration()=0; - virtual void AgressiveShelfSync()=0; + virtual void SoftMigration(double* sheet_ungrounding)=0; + virtual void ShelfSync()=0; + virtual void PotentialSheetUngrounding(Vec potential_sheet_ungrounding)=0; virtual void MigrateGroundingLine()=0; virtual int UpdateShelfStatus(Vec new_shelf_nodes)=0; virtual void UpdateShelfFlags(double* new_shelf_nodes)=0; + virtual int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf)=0; /*Implementation: */ Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 7322) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 7323) @@ -309,6 +309,11 @@ } /*}}}*/ -/*FUNCTION Penta::AgressiveShelfSync{{{1*/ -void Penta::AgressiveShelfSync(void){ +/*FUNCTION Penta::SoftMigration{{{1*/ +void Penta::SoftMigration(double* sheet_ungrounding){ + _error_("not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::ShelfSync{{{1*/ +void Penta::ShelfSync(void){ _error_("not supported yet!"); } @@ -5253,4 +5258,19 @@ } /*}}}*/ +/*FUNCTION Penta::IsNodeOnShelf {{{1*/ +bool Penta::IsNodeOnShelfFromFlags(double* flags){ + + int i; + bool shelf=false; + + for(i=0;i<6;i++){ + if (flags[nodes[i]->Sid()]){ + shelf=true; + break; + } + } + return shelf; +} +/*}}}*/ /*FUNCTION Penta::IsOnSurface{{{1*/ bool Penta::IsOnSurface(void){ @@ -5500,4 +5520,9 @@ *pnumvertices=NUMVERTICES; *pnumnodes=numnodes; +} +/*}}}*/ +/*FUNCTION Penta::PotentialSheetUngrounding{{{1*/ +void Penta::PotentialSheetUngrounding(Vec potential_sheet_ungrounding){ + _error_("not supported yet!"); } /*}}}*/ @@ -6150,2 +6175,7 @@ } /*}}}*/ +/*FUNCTION Penta::UpdatePotentialSheetUngrounding{{{1*/ +int Penta::UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf){ + _error_("Not implemented yet"); +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 7322) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 7323) @@ -110,5 +110,7 @@ void MaxVz(double* pmaxvz, bool process_units); void AgressiveMigration(); - void AgressiveShelfSync(); + void SoftMigration(double* sheet_ungrounding); + void PotentialSheetUngrounding(Vec potential_sheet_ungrounding); + void ShelfSync(); void MigrateGroundingLine(); void MinVel(double* pminvel, bool process_units); @@ -129,4 +131,5 @@ int UpdateShelfStatus(Vec new_shelf_nodes); void UpdateShelfFlags(double* new_shelf_nodes); + int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf); double TimeAdapt(); int* GetHorizontalNeighboorSids(void); @@ -239,4 +242,5 @@ bool IsOnShelf(void); bool IsNodeOnShelf(void); + bool IsNodeOnShelfFromFlags(double* flags); bool IsOnWater(void); double MinEdgeLength(double xyz_list[6][3]); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 7322) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 7323) @@ -335,55 +335,4 @@ } -/*}}}*/ -/*FUNCTION Tria::AgressiveShelfSync{{{1*/ -void Tria::AgressiveShelfSync(void){ - - - double *values = NULL; - double h[3],s[3],b[3],ba[3]; - double bed_hydro; - double rho_water,rho_ice,density; - int i; - bool elementonshelf = false; - - /*Recover info at the vertices: */ - Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); - if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); - - GetParameterListOnVertices(&h[0],ThicknessEnum); - GetParameterListOnVertices(&s[0],SurfaceEnum); - GetParameterListOnVertices(&b[0],BedEnum); - GetParameterListOnVertices(&ba[0],BathymetryEnum); - - /*material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - density=rho_ice/rho_water; - - /*go through vertices, and update inputs, considering them to be TriaVertex type: */ - for(i=0;i<3;i++){ - if(b[i]==ba[i]){ - - nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false)); - nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true)); - } - else{ - nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true)); - nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false)); - - } - } - - /*Now, update shelf status of element. An element can only be on shelf if all its nodes are on shelf: */ - elementonshelf=false; - for(i=0;i<3;i++){ - if(nodes[i]->IsOnShelf()){ - elementonshelf=true; - break; - } - } - this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf)); -} /*}}}*/ /*FUNCTION Tria::AverageOntoPartition {{{1*/ @@ -4348,4 +4297,19 @@ } /*}}}*/ +/*FUNCTION Tria::IsNodeOnShelf {{{1*/ +bool Tria::IsNodeOnShelfFromFlags(double* flags){ + + int i; + bool shelf=false; + + for(i=0;i<3;i++){ + if (flags[nodes[i]->Sid()]){ + shelf=true; + break; + } + } + return shelf; +} +/*}}}*/ /*FUNCTION Tria::IsOnWater {{{1*/ bool Tria::IsOnWater(){ @@ -4695,4 +4659,45 @@ } /*}}}*/ +/*FUNCTION Tria::PotentialSheetUngrounding{{{1*/ +void Tria::PotentialSheetUngrounding(Vec potential_sheet_ungrounding){ + + + double *values = NULL; + double h[3],s[3],b[3],ba[3]; + double bed_hydro; + double rho_water,rho_ice,density; + int i; + bool elementonshelf = false; + + /*Recover info at the vertices: */ + Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); + if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); + + GetParameterListOnVertices(&h[0],ThicknessEnum); + GetParameterListOnVertices(&s[0],SurfaceEnum); + GetParameterListOnVertices(&b[0],BedEnum); + GetParameterListOnVertices(&ba[0],BathymetryEnum); + + /*material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + density=rho_ice/rho_water; + + /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */ + for(i=0;i<3;i++){ + if (!nodes[i]->IsOnShelf()){ + + /*This node is on the sheet, near the grounding line. See if wants to unground. To + * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */ + bed_hydro=-density*h[i]; + if (bed_hydro>ba[i]){ + /*ok, this node wants to unground. flag it: */ + VecSetValue(potential_sheet_ungrounding,nodes[i]->Sid(),1,INSERT_VALUES); + } + } + } +} +/*}}}*/ /*FUNCTION Tria::ProcessResultsUnits{{{1*/ void Tria::ProcessResultsUnits(void){ @@ -4799,4 +4804,119 @@ } +/*}}}*/ +/*FUNCTION Tria::ShelfSync{{{1*/ +void Tria::ShelfSync(void){ + + + double *values = NULL; + double h[3],s[3],b[3],ba[3]; + double bed_hydro; + double rho_water,rho_ice,density; + int i; + bool elementonshelf = false; + + /*Recover info at the vertices: */ + Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); + if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); + + GetParameterListOnVertices(&h[0],ThicknessEnum); + GetParameterListOnVertices(&s[0],SurfaceEnum); + GetParameterListOnVertices(&b[0],BedEnum); + GetParameterListOnVertices(&ba[0],BathymetryEnum); + + /*material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + density=rho_ice/rho_water; + + /*go through vertices, and update inputs, considering them to be TriaVertex type: */ + for(i=0;i<3;i++){ + if(b[i]==ba[i]){ + + nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false)); + nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true)); + } + else{ + nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true)); + nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false)); + + } + } + + /*Now, update shelf status of element. An element can only be on shelf if all its nodes are on shelf: */ + elementonshelf=false; + for(i=0;i<3;i++){ + if(nodes[i]->IsOnShelf()){ + elementonshelf=true; + break; + } + } + this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf)); +} +/*}}}*/ +/*FUNCTION Tria::SoftMigration{{{1*/ +void Tria::SoftMigration(double* sheet_ungrounding){ + + + double *values = NULL; + double h[3],s[3],b[3],ba[3]; + double bed_hydro; + double rho_water,rho_ice,density; + int i; + bool elementonshelf = false; + + /*Recover info at the vertices: */ + Input* surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); + if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); + + GetParameterListOnVertices(&h[0],ThicknessEnum); + GetParameterListOnVertices(&s[0],SurfaceEnum); + GetParameterListOnVertices(&b[0],BedEnum); + GetParameterListOnVertices(&ba[0],BathymetryEnum); + + /*material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + density=rho_ice/rho_water; + + /*go through vertices, and update inputs, considering them to be TriaVertex type: */ + for(i=0;i<3;i++){ + if (nodes[i]->IsOnShelf()){ + /*This node is on the shelf. See if its bed is going under the bathymetry: */ + if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway. + /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */ + b[i]=ba[i]; + s[i]=b[i]+h[i]; + } + else{ + /*do nothing, we are still floating.*/ + } + } + else{ + /*This node is on the sheet, near the grounding line. See if wants to unground. To + * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */ + bed_hydro=-density*h[i]; + if (bed_hydro>ba[i]){ + + /*Now, are we connected to the grounding line, if so, go forward, otherwise, bail out: */ + if(sheet_ungrounding[nodes[i]->Sid()]){ + /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */ + s[i]=(1-density)*h[i]; + b[i]=-density*h[i]; + } + } + else{ + /*do nothing, we are still grounded.*/ + } + } + } + + /*Surface and bed are updated. Update inputs:*/ + surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i]; + bed_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=b[i]; + + } /*}}}*/ /*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ @@ -5548,2 +5668,20 @@ } /*}}}*/ +/*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{1*/ +int Tria::UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf){ + + /*intermediary: */ + int i; + int nflipped=0; + + /*Ok, go through our 3 nodes, and whoever is on the potential_sheet_ungrounding, ends up in nodes_on_iceshelf: */ + for(i=0;i<3;i++){ + if (potential_sheet_ungrounding[nodes[i]->Sid()]){ + VecSetValue(vec_nodes_on_iceshelf,nodes[i]->Sid(),1,INSERT_VALUES); + } + /*Figure out if we flipped: */ + if (potential_sheet_ungrounding[nodes[i]->Sid()] != nodes_on_iceshelf[nodes[i]->Sid()])nflipped++; + } + return nflipped; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 7322) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 7323) @@ -84,4 +84,5 @@ bool IsOnShelf(); bool IsNodeOnShelf(); + bool IsNodeOnShelfFromFlags(double* flags); bool IsOnWater(); void GetSolutionFromInputs(Vec solution); @@ -114,5 +115,7 @@ void MaxVz(double* pmaxvz, bool process_units); void AgressiveMigration(); - void AgressiveShelfSync(); + void SoftMigration(double* sheet_ungrounding); + void ShelfSync(); + void PotentialSheetUngrounding(Vec potential_sheet_ungrounding); void MigrateGroundingLine(); void MinVel(double* pminvel, bool process_units); @@ -132,4 +135,5 @@ void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); int UpdateShelfStatus(Vec new_shelf_nodes); + int UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf); void UpdateShelfFlags(double* new_shelf_nodes); double TimeAdapt();