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

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

/*Headers: {{{1*/
#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(){{{1*/
Elements::Elements(){
        return;
}
/*}}}*/
/*FUNCTION Elements::Elements(int in_enum){{{1*/
Elements::Elements(int in_enum): DataSet(in_enum){
        //do nothing;
        return;
}
/*}}}*/
/*FUNCTION Elements::~Elements(){{{1*/
Elements::~Elements(){
        return;
}
/*}}}*/



/*Object management*/
/*FUNCTION Elements::ProcessResultsUnits{{{1*/
void Elements::ProcessResultsUnits(void){

        int i;

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

        for (i=0;i<this->Size();i++){
                Element* element=(Element*)this->GetObjectByOffset(i);
                element->DeleteResults();
        }
}
/*}}}*/
/*FUNCTION Elements::ToResults{{{1*/
void Elements::ToResults(Results* results,Parameters* parameters,int step, double time){

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

        /*I/O strategy: */
        bool   io_gather=true; //the default

        /*Recover parameters: */
        parameters->FindParam(&io_gather,IoGatherEnum);


        /*create patch object out of all results in this dataset: */
        patch=this->ResultsToPatch();

        /*Gather onto master cpu 0, if needed: */
        #ifdef _PARALLEL_
        if(io_gather)patch->MPI_Gather();
        #endif
        
        /*create result object and add to results dataset:*/
        if (patch->maxvertices && patch->maxnodes){
                results->AddObject(new IntExternalResult(results->Size()+1,PatchVerticesEnum,patch->maxvertices,step,time));
                results->AddObject(new IntExternalResult(results->Size()+1,PatchNodesEnum,   patch->maxnodes,step,time));
                results->AddObject(new DoubleMatExternalResult(results->Size()+1,PatchEnum,patch->values,patch->numrows,patch->numcols,step,time));
        }

        /*Free ressources:*/
        delete patch;

}
/*}}}*/
/*FUNCTION Elements::ResultsToPatch{{{1*/
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=(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;
        }

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

        MPI_Reduce (&numnodes,&max_numnodes,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD );
        MPI_Bcast(&max_numnodes,1,MPI_INT,0,MPI_COMM_WORLD);
        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=(Element*)this->GetObjectByOffset(i);
                element->PatchFill(&count,patch);
        }

        return patch;

}
/*}}}*/