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

/*!\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 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      numberofvertices;
        int      nflipped,local_nflipped;
        double*  nodes_on_iceshelf             = NULL;
        double*  elements_touching_iceshelf    = NULL;
        Vec      vec_element_touching_iceshelf = NULL;
        Vec      vec_nodes_on_iceshelf         = NULL;
        Node*    node                          = NULL;
        Element* element                       = NULL;


        /*recover parameters: */
        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
        numberofvertices=vertices->NumberOfVertices();

        /*recover vec_nodes_on_iceshelf*/
        vec_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(analysis_type)){
                        if(node->IsFloating()){
                                VecSetValue(vec_nodes_on_iceshelf,node->Sid(),1.0,INSERT_VALUES);
                        }
                }
        }

        /*Assemble vector and serialize: */
        VecAssemblyBegin(vec_nodes_on_iceshelf);
        VecAssemblyEnd(vec_nodes_on_iceshelf);
        VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);

        nflipped=1; //bootstrap
        while(nflipped){
                
                /*Vector of size number of elements*/
                vec_element_touching_iceshelf=NewVec(elements->NumberOfElements(),true);

                /*Figure out if any of the nodes of the element will be floating -> elements neighbouting the floating ice*/
                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 and serialize: */
                VecAssemblyBegin(vec_element_touching_iceshelf);
                VecAssemblyEnd(vec_element_touching_iceshelf);
                VecToMPISerial(&elements_touching_iceshelf,vec_element_touching_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:*/
                VecFree(&vec_element_touching_iceshelf);
                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;
}
/*}}}*/