source: issm/trunk/src/c/classes/Patch.cpp@ 15396

/*!\file Patch.c
 * \brief: implementation of the Patch object
 */

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

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "./classes.h"
#include "../shared/shared.h"

/*Object constructors and destructors:*/
/*FUNCTION Patch::Patch() default constructor {{{*/
Patch::Patch(){
        this->numrows     = 0;
        this->numcols     = 0;
        this->maxvertices = 0;
        this->maxnodes    = 0;
        this->values      = NULL;
}
/*}}}*/
/*FUNCTION Patch::Patch(int numrows, int maxvertices, int maxnodes){{{*/
Patch::Patch(int in_numrows, int in_maxvertices, int in_maxnodes){

        this->numrows=in_numrows;
        this->maxvertices=in_maxvertices;
        this->maxnodes=in_maxnodes;
        this->numcols=1  //enum_type
                +1           //step
                +1           //time 
                +1           //element id
                +1           //interpolation type
                +maxvertices //vertices
                +maxnodes;   //nodes

        //Allocate values and fill with NaN:
        if (this->numcols*this->numrows==0){
                this->values=NULL;
        }
        else{
                this->values=xNew<IssmDouble>(this->numcols*this->numrows);
                for(int i=0;i<this->numrows;i++){
                        for(int j=0;j<this->numcols;j++){
                                this->values[i*this->numcols+j]=NAN;
                        }
                }
        }

}
/*}}}*/
/*FUNCTION Patch::~Patch(){{{*/
Patch::~Patch(){
        xDelete<IssmDouble>(values);
}
/*}}}*/

/*Object methods*/
/*FUNCTION Patch::fillelementinfo{{{*/
void Patch::fillelementinfo(int count, int element_id, int* vertices_ids, int num_vertices){

        int i;
        IssmDouble* row=NULL;

        /*point to the start of the row: */
        row=this->values+count*this->numcols;

        /*Let's remember what is on a row: 
         enum_type step time element_id interpolation vertices_ids nodal_values
         */
        row[3]=element_id;
        for(i=0;i<num_vertices;i++){
                row[5+i]=vertices_ids[i];
        }

}
/*}}}*/
/*FUNCTION Patch::fillresultinfo{{{*/
void Patch::fillresultinfo(int count,int enum_type,int step, IssmDouble time, int interpolation, IssmDouble* nodal_values, int num_nodes){

        int i;
        IssmDouble* row=NULL;

        /*point to the start of the row: */
        row=this->values+count*this->numcols;

        /*Let's remember what is on a row: 
         enum_type step time element_id interpolation vertices_ids nodal_values */
        row[0]=enum_type;
        row[1]=(IssmDouble)step;
        row[2]=time;
        row[4]=interpolation;
        for(i=0;i<num_nodes;i++){
                row[5+this->maxvertices+i]=nodal_values[i];
        }

}
/*}}}*/
/*FUNCTION Patch::Gather{{{*/
void Patch::Gather(void){

        int         count;
        int         my_rank;
        int         num_procs;
        int         total_numrows;
        int         node_numrows;
        IssmDouble     *total_values  = NULL;
        #ifdef _HAVE_MPI_
        MPI_Status  status;
        #endif

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

        /*First, figure out total number of rows combining all the cpus: */
        #ifdef _HAVE_MPI_
        MPI_Reduce(&this->numrows,&total_numrows,1,MPI_INT,MPI_SUM,0,IssmComm::GetComm() );
        MPI_Bcast(&total_numrows,1,MPI_INT,0,IssmComm::GetComm());
        #else
        total_numrows=this->numrows;
        #endif

        /*return if patch empty*/
        if(total_numrows==0) return;

        /*Now, allocate buffer to holds all the values, on node 0: */
        if(my_rank==0)total_values=xNew<IssmDouble>(this->numcols*total_numrows);

        /*Start by copying node 0 values onto total_values: */
        if(my_rank==0){
                count=0;
                xMemCpy<IssmDouble>(total_values+count,this->values,this->numcols*this->numrows);
                count+=this->numrows*this->numcols;
        }

        /*Now, ask other nodes to send their values: */
        #ifdef _HAVE_MPI_
        for(int i=1;i<num_procs;i++){
                if (my_rank==i){ 
                        MPI_Send(&this->numrows,1,MPI_INT,0,1,IssmComm::GetComm());   
                        if (this->numrows)MPI_Send(this->values,this->numrows*this->numcols,MPI_DOUBLE,0,1,IssmComm::GetComm()); 
                }
                if (my_rank==0){
                        MPI_Recv(&node_numrows,1,MPI_INT,i,1,IssmComm::GetComm(),&status); 
                        if (node_numrows)MPI_Recv(total_values+count,node_numrows*this->numcols,MPI_DOUBLE,i,1,IssmComm::GetComm(),&status);
                        count+=node_numrows*this->numcols;
                }
        }       
        #endif

        /*Now, node 0 has total_values, of size total_numrows*this->numcols. Update the fields in the patch, to reflect this new 
         * reality. For other cpus, no point in keeping their data anymore: */
        if(my_rank==0){
                this->numrows=total_numrows;
                xDelete<IssmDouble>(this->values);
                this->values=total_values;
        }
        #ifdef _HAVE_MPI_
        else{
                this->numrows=0;
                xDelete<IssmDouble>(this->values);
        }
        #endif
}/*}}}*/