source: issm/trunk/src/c/Container/Elements.cpp@ 13975

/*
 * \file Elements.c
 * \brief: implementation of the Elements class, derived from DataSet class
 */

/*Headers: {{{*/
#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <vector>
#include <functional>
#include <algorithm>
#include <iostream>

#include "./DataSet.h"
#include "../shared/shared.h"
#include "../include/include.h"
#include "../EnumDefinitions/EnumDefinitions.h"

using namespace std;
/*}}}*/

/*Object constructors and destructor*/
/*FUNCTION Elements::Elements(){{{*/
Elements::Elements(){
        enum_type=MeshElementsEnum;
        return;
}
/*}}}*/
/*FUNCTION Elements::~Elements(){{{*/
Elements::~Elements(){
        return;
}
/*}}}*/

/*Object management*/
/*FUNCTION Elements::Configure{{{*/
void Elements::Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){

        vector<Object*>::iterator object;
        Element* element=NULL;

        for ( object=objects.begin() ; object < objects.end(); object++ ){

                element=dynamic_cast<Element*>((*object));
                element->Configure(elements,loads,nodes,materials,parameters);

        }

}
/*}}}*/
/*FUNCTION Elements::ProcessResultsUnits{{{*/
void Elements::ProcessResultsUnits(void){

        //Process results to be output in the correct units
        for(int i=0;i<this->Size();i++){
                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
                element->ProcessResultsUnits();
        }
}
/*}}}*/
/*FUNCTION Elements::DeleteResults{{{*/
void Elements::DeleteResults(void){

        for (int i=0;i<this->Size();i++){
                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
                element->DeleteResults();
        }
}
/*}}}*/
/*FUNCTION Elements::ResultsToPatch{{{*/
Patch* Elements::ResultsToPatch(void){ 

        /*output: */
        Patch* patch=NULL;

        /*intermediary: */
        int i;
        int count;
        int numrows;
        int numvertices;
        int numnodes;
        int max_numvertices;
        int max_numnodes;
        int element_numvertices;
        int element_numrows;
        int element_numnodes;

        /*We are going to extract from the results within the elements, the desired results , and create a table 
         * of patch information, that will hold, for each element that computed the result that 
         * we desire, the enum_type of the result, the step and time, the id of the element, the interpolation type, the vertices ids, and the values 
         * at the nodes (could be different from the vertices). This will be used for visualization purposes. 
         * For example, we could build the following patch table, for velocities: 
         * 
         1. on a Beam element, Vx, at step 1, time .5, element id 1, interpolation type P0 (constant), vertices ids 1 and 2, one constant value 4.5
         VxEnum 1  .5  1 P0  1 2       4.5 NaN  NaN
         2. on a Tria element, Vz, at step 2, time .8, element id 2, interpolation type P1 (linear), vertices ids 1 3 and 4, with values at 3 nodes 4.5, 3.2, 2.5
         VzEnum 2  .8  2 P1  1 3 4     4.5 3.2  2.5
         * ... etc ...
         *
         * So what do we need to build the table: the maximum number of vertices included in the table, 
         * and the maximum number of nodal values, as well as the number of rows. Once we have that, 
         * we ask the elements to fill their own row in the table, by looping on the elememnts. 
         *
         * We will use the Patch object, which will store all of the information needed, and will be able 
         * to output itself to disk on its own. See object/Patch.h for format of this object.*/

        /*First, determine maximum number of vertices, nodes, and number of results: */
        numrows=0;
        numvertices=0;
        numnodes=0;

        for(i=0;i<this->Size();i++){

                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
                element->PatchSize(&element_numrows,&element_numvertices,&element_numnodes);

                numrows+=element_numrows;
                if(element_numvertices>numvertices)numvertices=element_numvertices;
                if(element_numnodes>numnodes)numnodes=element_numnodes;
        }

        /*Synchronize across cluster, so as to not end up with different sizes for each patch on each cpu: */
        #ifdef _HAVE_MPI_
        MPI_Reduce (&numvertices,&max_numvertices,1,MPI_INT,MPI_MAX,0,IssmComm::GetComm() );
        MPI_Bcast(&max_numvertices,1,MPI_INT,0,IssmComm::GetComm());
        numvertices=max_numvertices;

        MPI_Reduce (&numnodes,&max_numnodes,1,MPI_INT,MPI_MAX,0,IssmComm::GetComm() );
        MPI_Bcast(&max_numnodes,1,MPI_INT,0,IssmComm::GetComm());
        numnodes=max_numnodes;
        #endif

        /*Ok, initialize Patch object: */
        patch=new Patch(numrows,numvertices,numnodes);

        /*Now, go through elements, and fill the Patch object: */
        count=0;
        for(i=0;i<this->Size();i++){
                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
                element->PatchFill(&count,patch);
        }

        return patch;
}
/*}}}*/
/*FUNCTION Elements::SetCurrentConfiguration{{{*/
void Elements::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){

        vector<Object*>::iterator object;
        Element* element=NULL;

        for ( object=objects.begin() ; object < objects.end(); object++ ){

                element=dynamic_cast<Element*>((*object));
                element->SetCurrentConfiguration(elements,loads,nodes,materials,parameters);

        }

}
/*}}}*/
/*FUNCTION Elements::ToResults{{{*/
void Elements::ToResults(Results* results,Parameters* parameters){

        int my_rank;
        int num_procs;

        Patch               *patch              = NULL;
        int                 *resultsenums       = NULL;
        int                 *resultssizes       = NULL;
        int                 *resultssteps       = NULL;
        IssmDouble          *resultstimes       = NULL;
        IssmDouble          *vector_serial      = NULL;
        Vector<IssmDouble> *vector = NULL;
        bool                io_gather;
        bool                results_as_patches;
        int                 numberofvertices     ,numberofelements;
        int                 numberofresults      ,vectorsize;
        int                 rank;
        int                 minrank;

        /*recover my_rank:*/
        my_rank=IssmComm::GetRank();
        num_procs=IssmComm::GetSize();

        /*Recover parameters: */
        parameters->FindParam(&io_gather,SettingsIoGatherEnum);
        parameters->FindParam(&results_as_patches,SettingsResultsAsPatchesEnum);
        parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
        parameters->FindParam(&numberofelements,MeshNumberofelementsEnum);

        if(!results_as_patches){
                /*No patch here, we prepare vectors*/

                /*Get rank of first cpu that has results*/
                #ifdef _HAVE_MPI_
                if(this->Size()) rank=my_rank;
                else rank=num_procs;
                MPI_Allreduce (&rank,&minrank,1,MPI_INT,MPI_MIN,IssmComm::GetComm());
                #else
                minrank=my_rank;
                #endif

                /*see what the first element of this partition has in stock (this is common to all partitions)*/
                if(my_rank==minrank){
                        if(this->Size()==0) _error_("Cannot write results because there is no element??");
                        Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(0));
                        element->ListResultsInfo(&resultsenums,&resultssizes,&resultstimes,&resultssteps,&numberofresults);
                }
                #ifdef _HAVE_MPI_
                MPI_Bcast(&numberofresults,1,MPI_DOUBLE,minrank,IssmComm::GetComm());
                #endif

                /*Get out if there is no results. Otherwise broadcast info*/
                if(!numberofresults) return;
                #ifdef _HAVE_MPI_
                if(my_rank!=minrank){
                        resultsenums=xNew<int>(numberofresults);
                        resultssizes=xNew<int>(numberofresults);
                        resultstimes=xNew<IssmDouble>(numberofresults);
                        resultssteps=xNew<int>(numberofresults);
                }
                MPI_Bcast(resultsenums,numberofresults,MPI_INT,minrank,IssmComm::GetComm());
                MPI_Bcast(resultssizes,numberofresults,MPI_INT,minrank,IssmComm::GetComm());
                MPI_Bcast(resultstimes,numberofresults,MPI_DOUBLE,minrank,IssmComm::GetComm());
                MPI_Bcast(resultssteps,numberofresults,MPI_INT,minrank,IssmComm::GetComm());
                #endif

                /*Loop over all results and get nodal vector*/
                for(int i=0;i<numberofresults;i++){

                        /*Get vector for result number i*/
                        if(resultssizes[i]==P1Enum)      vectorsize=numberofvertices;
                        else if(resultssizes[i]==P0Enum) vectorsize=numberofelements;
                        else _error_("Unkown result size: " << EnumToStringx(resultssizes[i]));
                        vector=new Vector<IssmDouble>(vectorsize);

                        for(int j=0;j<this->Size();j++){
                                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(j));
                                element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]);
                        }
                        vector->Assemble();

                        /*Serialize and add to results*/
                        vector_serial=vector->ToMPISerial();
                        if(my_rank==0){
                                /*No need to add this vector for all cpus*/
                                #ifdef _HAVE_ADOLC_
                                IssmPDouble* vector_serial_passive=xNew<IssmPDouble>(vectorsize);
                                for(int k=0;k<vectorsize;k++)vector_serial_passive[k]=reCast<IssmPDouble>(vector_serial[k]);
                                results->AddObject(new GenericExternalResult<double*>(results->Size()+1,resultsenums[i],vector_serial_passive,vectorsize,1,resultssteps[i],resultstimes[i]));
                                xDelete<IssmPDouble>(vector_serial_passive);
                                #else
                                results->AddObject(new GenericExternalResult<double*>(results->Size()+1,resultsenums[i],vector_serial,vectorsize,1,resultssteps[i],resultstimes[i]));
                                #endif
                        }

                        /*clean up*/
                        xdelete(&vector);
                        xDelete<IssmDouble>(vector_serial);
                }
        }
        else{
                /*create patch object out of all results in this dataset: */
                patch=this->ResultsToPatch();

                /*Gather onto master cpu 0, if needed: */
                if(io_gather)patch->Gather();

                /*create result object and add to results dataset:*/
                if (patch->maxvertices && patch->maxnodes){
                        results->AddObject(new GenericExternalResult<int>(results->Size()+1,PatchVerticesEnum,patch->maxvertices,1,0));
                        results->AddObject(new GenericExternalResult<int>(results->Size()+1,PatchNodesEnum,   patch->maxnodes,1,0));
                        #ifdef _HAVE_ADOLC_
                        IssmPDouble* values_passive=xNew<IssmPDouble>(patch->numrows*patch->numcols);
                        for(int k=0;k<(patch->numrows*patch->numcols);k++)values_passive[k]=reCast<IssmPDouble>(patch->values[k]);
                        results->AddObject(new GenericExternalResult<double*>(results->Size()+1,PatchEnum,  values_passive,patch->numrows,patch->numcols,1,0));
                        xDelete<IssmPDouble>(values_passive);
                        #else
                        results->AddObject(new GenericExternalResult<double*>(results->Size()+1,PatchEnum,  patch->values,patch->numrows,patch->numcols,1,0));
                        #endif
                }
        }

        /*Free ressources:*/
        xDelete<int>(resultsenums);
        xDelete<int>(resultssizes);
        xDelete<int>(resultssteps);
        xDelete<IssmDouble>(resultstimes);
        delete patch;
}
/*}}}*/
/*FUNCTION Elements::MaxNumNodes{{{*/
int Elements::MaxNumNodes(void){

        int max=0;
        int allmax;
        int numnodes=0;

        /*Now go through all elements, and get how many nodes they own, unless they are clone nodes: */
        for(int i=0;i<this->Size();i++){

                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
                numnodes=element->GetNumberOfNodes();
                if(numnodes>max)max=numnodes;
        }

        /*Grab max of all cpus: */
        #ifdef _HAVE_MPI_
        MPI_Allreduce((void*)&max,(void*)&allmax,1,MPI_INT,MPI_MAX,IssmComm::GetComm());
        max=allmax;
        #endif

        return max;
}
/*}}}*/
/*FUNCTION Elements::NumberOfElements{{{*/
int Elements::NumberOfElements(void){

        int local_nelem;
        int numberofelements;

        local_nelem=this->Size();
        #ifdef _HAVE_MPI_
        MPI_Allreduce ( (void*)&local_nelem,(void*)&numberofelements,1,MPI_INT,MPI_SUM,IssmComm::GetComm());
        #else
        numberofelements=local_nelem;
        #endif

        return numberofelements;
}
/*}}}*/
/*FUNCTION Elements::InputCopy{{{*/
void Elements::InputDuplicate(int input_enum,int output_enum){

        for(int i=0;i<this->Size();i++){
                Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
                element->InputDuplicate(input_enum,output_enum);
        }
}
/*}}}*/