Index: /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.cpp (revision 10292) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.cpp (revision 10292) @@ -0,0 +1,30 @@ +/*!\file GroundingLineMigrationx + * \brief: migration grounding line position. + */ + +#include "./GroundingLineMigrationx.h" +#include "./GroundingLineMigrationxLocal.h" + +#include "../../shared/shared.h" +#include "../../io/io.h" +#include "../../include/include.h" +#include "../../io/io.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../Container/Container.h" + +void GroundingLineMigrationx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ + + int migration_style; + + _printf_(VerboseModule()," Migrating grounding line\n"); + + /*retrieve parameters: */ + parameters->FindParam(&migration_style,GroundinglineMigrationEnum); + + /*call different migration modules, according to user wishes: */ + if(migration_style==AgressiveMigrationEnum) 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_(true,"%s\n","NoneEnum supplied for migration style, doing nothing!"); + else _error_("%s not supported yet!",EnumToStringx(migration_style)); +} Index: /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.h =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.h (revision 10292) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.h (revision 10292) @@ -0,0 +1,13 @@ +/*!\file: GroundingLineMigrationx.h + * \brief header file for Grounding Line Migration + */ + +#ifndef _GROUNDINGLINEMIGRATIONX_H +#define _GROUNDINGLINEMIGRATIONX_H + +#include "../../Container/Container.h" + +/* local prototypes: */ +void GroundingLineMigrationx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); + +#endif /* _GROUNDINGLINEMIGRATIONX_H */ Index: /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxLocal.h =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxLocal.h (revision 10292) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxLocal.h (revision 10292) @@ -0,0 +1,22 @@ +/*!\file: GroundingLineMigrationxLocal.h + * \brief local header file for Grounding Line Migration + */ + +#ifndef _GROUNDINGLINEMIGRATIONXLOCAL_H +#define _GROUNDINGLINEMIGRATIONXLOCAL_H + +class Elements; +class Nodes; +class Parameters; + +/* 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* 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 10292) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxUtils.cpp (revision 10292) @@ -0,0 +1,353 @@ +/*!\file GroundingLineMigrationx + * \brief: migration grounding line position. + */ + +#include "./GroundingLineMigrationx.h" +#include "./GroundingLineMigrationxLocal.h" +#include "../../shared/shared.h" +#include "../../io/io.h" +#include "../../io/io.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +/*FUNCTION AgressiveMigration{{{1*/ +void AgressiveMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ + /*Here, whatever nodes inside the ice sheet want to unground, we allow -> instantaneous transmission of water through the bedrock: */ + + int i; + Element* element=NULL; + + _printf_(VerboseModule()," Migrating grounding line\n"); + + /*Carry out grounding line migration for those elements: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->AgressiveMigration(); + } + + /*Synchronize shelf status: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->ShelfSync(); + } + +} +/*}}}*/ +/*FUNCTION SoftMigration {{{1*/ +void SoftMigration(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ + + /*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); + +} +/*}}}*/ +/*FUNCTION PotentialSheetUngrounding {{{1*/ +double* PotentialSheetUngrounding(Elements* elements,Nodes* nodes,Parameters* parameters){ + + int i; + Element* element=NULL; + Vec vec_potential_sheet_ungrounding=NULL; + double* potential_sheet_ungrounding=NULL; + int numnods; + int configuration_type; + + /*recover parameters: */ + parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */ + numnods=nodes->NumberOfNodes(configuration_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; +} +/*}}}*/ +/*FUNCTION PropagateShelfIntoConnexIceSheet {{{1*/ +double* PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding){ + + 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()]){ + local_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); + + /*Avoid leaks: */ + xfree((void**)&nodes_on_iceshelf); + xfree((void**)&elements_touching_iceshelf); + + /*Assemble:*/ + VecAssemblyBegin(vec_nodes_on_iceshelf); + VecAssemblyEnd(vec_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; +} +/*}}}*/ +/*FUNCTION CreateElementOnGroundingLine {{{1*/ +bool* CreateElementOnGroundingLine(Elements* elements,double* element_on_iceshelf){ + + int i; + int j; + Element *element = NULL; + bool *element_on_gl = NULL; + bool ongl=false; + int *neighboorsids = NULL; + int sid; + int shelf; + + /*Go through elements, and look for elements that can possibly have a grounding line migration. These + * are elements that touch the grounding line: */ + element_on_gl=(bool*)xmalloc(elements->Size()*sizeof(bool)); + + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + + sid=element->Sid(); + neighboorsids=element->GetHorizontalNeighboorSids(); + + shelf=(int)element_on_iceshelf[sid]; + ongl=false; + for(j=0;j<3;j++){ + if (neighboorsids[j]<0)continue; //this neighboor does not exist + if ((shelf==1) & (element_on_iceshelf[neighboorsids[j]]==1))continue; //both this element and this neighboor are on th ice shelf + if ((shelf==0) & (element_on_iceshelf[neighboorsids[j]]==0))continue; //both this element and this neighboor are on the ice sheet + ongl=true; //neighboor j is on a different location than us, ie we are touching the grounding line. + } + element_on_gl[i]=ongl; + } + + return element_on_gl; +} +/*}}}*/ +/*FUNCTION CreateElementOnIceShelf {{{1*/ +double* CreateElementOnIceShelf(Elements* elements){ + + int i; + Element* element=NULL; + Vec vec_element_on_iceshelf=NULL; + double* element_on_iceshelf=NULL; + + /*Create vector holding all the elements IsFloating flags: */ + vec_element_on_iceshelf=NewVec(elements->NumberOfElements(),true); + + /*Loop through elements, and fill vec_element_on_iceshelf: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + VecSetValue(vec_element_on_iceshelf,element->Sid(),(int)element->IsFloating(),INSERT_VALUES); + } + + /*Assemble vector: */ + VecAssemblyBegin(vec_element_on_iceshelf); + VecAssemblyEnd(vec_element_on_iceshelf); + + /*Serialize vector: */ + VecToMPISerial(&element_on_iceshelf,vec_element_on_iceshelf); + + /*free ressouces: */ + VecFree(&vec_element_on_iceshelf); + + return element_on_iceshelf; +} +/*}}}*/ +/*FUNCTION CreateElementTouchingIceShelf {{{1*/ +double* CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf){ + + int i; + Element* element=NULL; + Vec vec_element_touching_iceshelf=NULL; + double* nodes_on_iceshelf=NULL; + + /*output: */ + double* element_touching_iceshelf=NULL; + + /*serialize: */ + VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf); + + /*Create vector holding all the elements IsFloating flags: */ + vec_element_touching_iceshelf=NewVec(elements->NumberOfElements(),true); + + /*Loop through elements, and fill vec_element_touching_iceshelf: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_iceshelf)?1.0:0.0,INSERT_VALUES); + } + + /*Assemble vector: */ + VecAssemblyBegin(vec_element_touching_iceshelf); + VecAssemblyEnd(vec_element_touching_iceshelf); + + /*Serialize vector: */ + VecToMPISerial(&element_touching_iceshelf,vec_element_touching_iceshelf); + + /*free ressouces: */ + VecFree(&vec_element_touching_iceshelf); + xfree((void**)&nodes_on_iceshelf); + + return element_touching_iceshelf; +} +/*}}}*/ +/*FUNCTION UpdateShelfStatus {{{1*/ +int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){ + + int i; + Element* element=NULL; + int configuration_type; + int numnods; + Vec vec_new_shelf_nodes=NULL; + double* new_shelf_nodes=NULL; + + /*output: */ + int local_nflipped=0; + int nflipped=0; + + /*recover parameters: */ + parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + /*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */ + numnods=nodes->NumberOfNodes(configuration_type); + vec_new_shelf_nodes=NewVec(numnods); + + /*Ok, now go through the elements that have gone through grounding line migration, + * and update their flags: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + if(element_touching_iceshelf[element->Sid()]){ + local_nflipped+=element->UpdateShelfStatus(vec_new_shelf_nodes); + } + } + MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); + + /*Serialize vec_new_shelf_nodes: */ + VecToMPISerial(&new_shelf_nodes,vec_new_shelf_nodes); + + /*Now, go through ALL elements, and update the status of the nodes, to propagate what happened at the grounding line: */ + /*Carry out grounding line migration for those elements: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->UpdateShelfFlags(new_shelf_nodes); + } + + /*Free ressources: */ + VecFree(&vec_new_shelf_nodes); + xfree((void**)&new_shelf_nodes); + + return nflipped; +} +/*}}}*/ +/*FUNCTION CreateNodesOnIceShelf {{{1*/ +Vec CreateNodesOnIceShelf(Nodes* nodes,int configuration_type){ + + /*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(configuration_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->InAnalysis(configuration_type)){ + if(node->IsFloating()){ + VecSetValue(nodes_on_iceshelf,node->Sid(),1.0,INSERT_VALUES); + } + } + } + + /*Assemble vector: */ + VecAssemblyBegin(nodes_on_iceshelf); + VecAssemblyEnd(nodes_on_iceshelf); + + return nodes_on_iceshelf; +} +/*}}}*/ Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 10291) +++ /issm/trunk/src/c/modules/modules.h (revision 10292) @@ -32,5 +32,5 @@ #include "./GetVectorFromInputsx/GetVectorFromInputsx.h" #include "./Gradjx/Gradjx.h" -#include "./GroundingLineMigrationx/GroundingLineMigrationx.h" +#include "./GroundinglineMigrationx/GroundinglineMigrationx.h" #include "./HoleFillerx/HoleFillerx.h" #include "./IceVolumex/IceVolumex.h"