Index: /issm/trunk/src/c/Container/DataSet.cpp =================================================================== --- /issm/trunk/src/c/Container/DataSet.cpp (revision 10299) +++ /issm/trunk/src/c/Container/DataSet.cpp (revision 10300) @@ -388,4 +388,12 @@ dataset->AddObject(doubleexternalresult);} break; + #ifdef _HAVE_GROUNDINGLINE_ + case BoolElementResultEnum:{ + BoolElementResult* boolelementresult=NULL; + boolelementresult=new BoolElementResult(); + boolelementresult->Demarshall(&marshalled_dataset); + dataset->AddObject(boolelementresult);} + break; + #endif default: _error_("could not recognize enum type: %s",EnumToStringx(enum_type)); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 10299) +++ /issm/trunk/src/c/Makefile.am (revision 10300) @@ -535,9 +535,9 @@ ./solutions/bedslope_core.cpp #}}} -#GroundingLine sources {{{1 -groundingline_sources= ./modules/GroundingLineMigrationx/GroundingLineMigrationx.cpp\ - ./modules/GroundingLineMigrationx/GroundingLineMigrationx.h\ - ./modules/GroundingLineMigrationx/GroundingLineMigrationxLocal.h\ - ./modules/GroundingLineMigrationx/GroundingLineMigrationxUtils.cpp +#Groundingline sources {{{1 +groundingline_sources= ./modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp\ + ./modules/GroundinglineMigrationx/GroundinglineMigrationx.h\ + ./modules/GroundinglineMigrationx/GroundinglineMigrationxLocal.h\ + ./modules/GroundinglineMigrationx/GroundinglineMigrationxUtils.cpp #}}} #Rifts sources {{{1 Index: sm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.cpp (revision 10299) +++ (revision ) @@ -1,30 +1,0 @@ -/*!\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: sm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.h =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationx.h (revision 10299) +++ (revision ) @@ -1,13 +1,0 @@ -/*!\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: sm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxLocal.h =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxLocal.h (revision 10299) +++ (revision ) @@ -1,22 +1,0 @@ -/*!\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: sm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxUtils.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundingLineMigrationxUtils.cpp (revision 10299) +++ (revision ) @@ -1,353 +1,0 @@ -/*!\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/GroundinglineMigrationx/GroundinglineMigrationx.cpp =================================================================== --- /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 10300) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 10300) @@ -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 10300) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.h (revision 10300) @@ -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 10300) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationxLocal.h (revision 10300) @@ -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 10300) +++ /issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationxUtils.cpp (revision 10300) @@ -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 CreateElementOnGroundingnine {{{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/solutions/transient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core.cpp (revision 10299) +++ /issm/trunk/src/c/solutions/transient_core.cpp (revision 10300) @@ -24,5 +24,5 @@ /*parameters: */ double finaltime,dt,yts; - bool control_analysis,isdiagnostic,isprognostic,isthermal; + bool control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline; bool time_adapt=false; int solution_type; @@ -45,8 +45,9 @@ femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum); femmodel->parameters->FindParam(&time_adapt,TimesteppingTimeAdaptEnum); - femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); femmodel->parameters->FindParam(&isdiagnostic,TransientIsdiagnosticEnum); femmodel->parameters->FindParam(&isprognostic,TransientIsprognosticEnum); femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); + femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); + if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum); @@ -86,5 +87,4 @@ _error_("ISSM was not compiled with diagnostic capabilities. Exiting"); #endif - } @@ -96,13 +96,12 @@ } - if (groundingline_migration!=NoneEnum){ + if(isgroundingline){ if(dim==3) _error_("Grounding line migration not implemented in 3d"); _printf_(VerboseSolution()," computing new grounding line position\n"); - #ifdef _HAVE_DIAGNOSTIC_ - GroundingLineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + #ifdef _HAVE_GROUNDINGLINE_ + GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); #else _error_("ISSM was not compiled with grounding line migration capabilities. Exiting"); #endif - }