source: issm/trunk/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationxUtils.cpp@ 10364

/*!\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 "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"

/*FUNCTION CreateElementOnGroundingline {{{1*/
bool*      CreateElementOnGroundingline(Elements* elements,double* element_on_iceshelf){ 

        bool    *element_on_gl = NULL;
        bool     ongl=false;
        int      i,j,sid,shelf;
        int     *neighboorsids = NULL;
        Element *element       = NULL;

        /*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;i<elements->Size();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;
        double  *element_on_iceshelf     = NULL;
        Element *element                 = NULL;
        Vec      vec_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;i<elements->Size();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,double* nodes_on_iceshelf){ 

        int      i;
        double  *element_touching_iceshelf     = NULL;
        Element *element                       = NULL;
        Vec      vec_element_touching_iceshelf = NULL;

        /*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;i<elements->Size();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);

        return element_touching_iceshelf;
}
/*}}}*/
/*FUNCTION CreateNodesOnIceShelf {{{1*/
Vec CreateNodesOnIceShelf(Nodes* nodes,Vertices* vertices,int configuration_type){ 

        int     i,numberofvertices;
        Vec     nodes_on_iceshelf = NULL;
        Node*   node              = NULL;

        /*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */
        numberofvertices=vertices->NumberOfVertices();
        nodes_on_iceshelf=NewVec(numberofvertices);

        /*Loop through nodes, and fill nodes_on_iceshelf: */
        for(i=0;i<nodes->Size();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;
}
/*}}}*/
/*FUNCTION PotentialSheetUngrounding {{{1*/
double*    PotentialSheetUngrounding(Elements* elements,Vertices* vertices,Parameters* parameters){ 

        int      i,numberofvertices;
        double*  vertices_potentially_ungrounding      = NULL;
        Vec      vec_vertices_potentially_ungrounding  = NULL;
        Element* element                               = NULL;

        /*Initialize vector with number of vertices*/
        numberofvertices=vertices->NumberOfVertices();
        vec_vertices_potentially_ungrounding=NewVec(numberofvertices); //grounded vertex that could start floating

        /*Fill vector vertices_potentially_floating: */
        for(i=0;i<elements->Size();i++){
                element=(Element*)elements->GetObjectByOffset(i);
                element->PotentialSheetUngrounding(vec_vertices_potentially_ungrounding);
        }

        /*Assemble vector: */
        VecAssemblyBegin(vec_vertices_potentially_ungrounding);
        VecAssemblyEnd(vec_vertices_potentially_ungrounding);

        /*Serialize vector: */
        VecToMPISerial(&vertices_potentially_ungrounding,vec_vertices_potentially_ungrounding);

        /*free ressouces and return: */
        VecFree(&vec_vertices_potentially_ungrounding);
        return vertices_potentially_ungrounding;
}
/*}}}*/
/*FUNCTION PropagateShelfIntoConnexIceSheet {{{1*/
double*    PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Vertices* vertices,Parameters* parameters,double* vertices_potentially_ungrounding){ 

        int      i,analysis_type;
        int      nflipped,local_nflipped;
        double*  nodes_on_iceshelf          = NULL;
        double*  elements_touching_iceshelf = NULL;
        Vec      vec_nodes_on_iceshelf      = NULL;
        Element* element                    = NULL;

        /*recover parameters: */
        parameters->FindParam(&analysis_type,AnalysisTypeEnum);

        /*recover vec_nodes_on_iceshelf and serialize: */
        vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,vertices,analysis_type);
        VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);

        nflipped=1; //bootstrap
        while(nflipped){
                
                /*get a list of potential elements that have nodes on ice shelf: */
                elements_touching_iceshelf=CreateElementTouchingIceShelf(elements,nodes_on_iceshelf);

                /*Go through elements_touching_iceshelf, and flag their nodes that can unground inside potential_sheet_ungrounding*/
                local_nflipped=0;
                for(i=0;i<elements->Size();i++){
                        element=(Element*)elements->GetObjectByOffset(i);
                        if(elements_touching_iceshelf[element->Sid()]){
                                local_nflipped+=element->UpdatePotentialSheetUngrounding(vertices_potentially_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**)&elements_touching_iceshelf);
                xfree((void**)&nodes_on_iceshelf);

                /*Assemble and serialize:*/
                VecAssemblyBegin(vec_nodes_on_iceshelf);
                VecAssemblyEnd(vec_nodes_on_iceshelf);
                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 UpdateShelfStatus {{{1*/
int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){

        int      i,numnods,configuration_type;
        int      nflipped=0;
        int      local_nflipped=0;
        Vec      vec_new_shelf_nodes = NULL;
        double*  new_shelf_nodes     = NULL;
        Element* element             = NULL;

        /*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;i<elements->Size();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;i<elements->Size();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;
}
/*}}}*/