source: issm/trunk/src/c/objects/Inputs/TriaVertexForcing.cpp@ 8375

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

#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 "../objects.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../shared/shared.h"
#include "../../Container/Container.h"
#include "../../include/include.h"

/*TriaVertexForcing constructors and destructor*/
/*FUNCTION TriaVertexForcing::TriaVertexForcing(){{{1*/
TriaVertexForcing::TriaVertexForcing(){
        enum_type=NoneEnum;
        values=NULL;
        numtimesteps=0;
        parameters=NULL;
        return;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* time0values,double time0,int numtimesteps,Parameters* parameters){{{1*/
TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* time0values,double time0,int in_numtimesteps,Parameters* parameters)
        :TriaRef(1)
{

        /*Set TriaRef*/
        this->SetElementType(P1Enum,0);
        this->element_type=P1Enum;

        /*Set Enum*/
        enum_type=in_enum_type;

        /*Allocate values and timesteps, and set first time step values: */
        this->numtimesteps=in_numtimesteps;
        this->values=(double*)xmalloc(this->numtimesteps*3*sizeof(double));
        this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
        this->values[0]=time0values[0];
        this->values[1]=time0values[1];
        this->values[2]=time0values[2];
        this->timesteps[0]=time0;

        this->parameters=parameters;

}
/*}}}*/
/*FUNCTION TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* in_values,double* in_time,int numtimesteps,Parameters* parameters){{{1*/
TriaVertexForcing::TriaVertexForcing(int in_enum_type,double* in_values,double* in_time,int in_numtimesteps,Parameters* parameters)
        :TriaRef(1)
{

        /*Set TriaRef*/
        this->SetElementType(P1Enum,0);
        this->element_type=P1Enum;

        /*Set Enum*/
        enum_type=in_enum_type;

        /*Allocate values and timesteps, and copy: */
        this->numtimesteps=in_numtimesteps;
        this->values=(double*)xmalloc(this->numtimesteps*3*sizeof(double));
        this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
        
        memcpy(this->values,in_values,numtimesteps*3*sizeof(double));
        memcpy(this->timesteps,in_time,numtimesteps*sizeof(double));

        this->parameters=parameters;

}
/*}}}*/
/*FUNCTION TriaVertexForcing::~TriaVertexForcing(){{{1*/
TriaVertexForcing::~TriaVertexForcing(){
        xfree((void**)&this->values);
        xfree((void**)&this->timesteps);
        parameters=NULL;
        return;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::AddTimeValues(double* values,int step,double time);{{{1*/
void TriaVertexForcing::AddTimeValues(double* values,int step,double time){

        /*insert values at time step: */
        if (step<0)_error_("timestep should start at least at time 0!");
        if (step>numtimesteps-1)_error_("timestep cannot be more than number of time steps");

        /*go ahead and plug: */
        this->timesteps[step]=time;
        this->values[3*step+0]=values[0];
        this->values[3*step+1]=values[1];
        this->values[3*step+2]=values[2];
}
/*}}}*/

/*Object virtual functions definitions:*/
/*FUNCTION TriaVertexForcing::Echo {{{1*/
void TriaVertexForcing::Echo(void){
        this->DeepEcho();
}
/*}}}*/
/*FUNCTION TriaVertexForcing::DeepEcho{{{1*/
void TriaVertexForcing::DeepEcho(void){

        int i;
        
        printf("TriaVertexForcing:\n");
        printf("   enum: %i (%s)\n",this->enum_type,EnumToStringx(this->enum_type));
        printf("   numtimesteps: %i\n",this->numtimesteps);
        for(i=0;i<this->numtimesteps;i++){
                printf("   time: %g  values: [%g %g %g]\n",this->timesteps[i],this->values[3*i+0],this->values[3*i+1],this->values[3*i+2]);
        }
}
/*}}}*/
/*FUNCTION TriaVertexForcing::Id{{{1*/
int    TriaVertexForcing::Id(void){ return -1; }
/*}}}*/
/*FUNCTION TriaVertexForcing::MyRank{{{1*/
int    TriaVertexForcing::MyRank(void){ 
        extern int my_rank;
        return my_rank; 
}
/*}}}*/
/*FUNCTION TriaVertexForcing::Marshall{{{1*/
void  TriaVertexForcing::Marshall(char** pmarshalled_dataset){

        char* marshalled_dataset=NULL;
        int   enum_value=0;

        /*recover marshalled_dataset: */
        marshalled_dataset=*pmarshalled_dataset;

        /*get enum value of TriaVertexForcing: */
        enum_value=TriaVertexForcingEnum;
        
        /*marshall enum: */
        memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
        
        /*marshall TriaVertexForcing data: */
        memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
        memcpy(marshalled_dataset,&numtimesteps,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps);
        memcpy(marshalled_dataset,values,numtimesteps*3*sizeof(double));marshalled_dataset+=numtimesteps*3*sizeof(double);
        memcpy(marshalled_dataset,timesteps,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double);

        *pmarshalled_dataset=marshalled_dataset;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::MarshallSize{{{1*/
int   TriaVertexForcing::MarshallSize(){
        
        return 
                +sizeof(enum_type)+
                +sizeof(numtimesteps)+
                3*numtimesteps*sizeof(double)+
                numtimesteps*sizeof(double)+
                +sizeof(int); //sizeof(int) for enum value
}
/*}}}*/
/*FUNCTION TriaVertexForcing::Demarshall{{{1*/
void  TriaVertexForcing::Demarshall(char** pmarshalled_dataset){

        char* marshalled_dataset=NULL;
        int   i;

        /*recover marshalled_dataset: */
        marshalled_dataset=*pmarshalled_dataset;

        /*this time, no need to get enum type, the pointer directly points to the beginning of the 
         *object data (thanks to DataSet::Demarshall):*/
        memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
        memcpy(&numtimesteps,marshalled_dataset,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps);
        /*allocate: */
        timesteps=(double*)xmalloc(numtimesteps*sizeof(double));
        values=(double*)xmalloc(numtimesteps*3*sizeof(double));
        
        memcpy(values,marshalled_dataset,3*numtimesteps*sizeof(double));marshalled_dataset+=3*numtimesteps*sizeof(double);
        memcpy(timesteps,marshalled_dataset,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double);

        /*return: */
        *pmarshalled_dataset=marshalled_dataset;
        return;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::Enum{{{1*/
int TriaVertexForcing::Enum(void){

        return TriaVertexForcingEnum;

}
/*}}}*/
/*FUNCTION TriaVertexForcing::copy{{{1*/
Object* TriaVertexForcing::copy() {
        
        return new TriaVertexForcing(this->enum_type,this->values,this->timesteps,this->numtimesteps,this->parameters);

}
/*}}}*/
        
/*TriaVertexForcing management*/
/*FUNCTION TriaVertexForcing::EnumType{{{1*/
int TriaVertexForcing::EnumType(void){

        return this->enum_type;

}
/*}}}*/
/*FUNCTION TriaVertexForcing::SpawnTriaInput{{{1*/
Input* TriaVertexForcing::SpawnTriaInput(int* indices){

        /*output*/
        TriaVertexForcing* outforcing=NULL;

        /*Create new Tria forcing (copy of current forcing)*/
        outforcing=new TriaVertexForcing(this->enum_type,this->values,this->timesteps,this->numtimesteps,this->parameters);

        /*Assign output*/
        return outforcing;

}
/*}}}*/
/*FUNCTION TriaVertexForcing::SpawnResult{{{1*/
ElementResult* TriaVertexForcing::SpawnResult(int step, double time){

        double triavalues[3];

        /*Ok, we want to spawn a TriaVertexElementResult. We have the time, just get 
         *the correct values at the three nodes: */
        this->GetTimeValues(&triavalues[0],time);

        return new TriaVertexElementResult(this->enum_type,&triavalues[0],step,time);

}
/*}}}*/

/*Object functions*/
/*FUNCTION TriaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
void TriaVertexForcing::GetParameterValue(double* pvalue,GaussTria* gauss){
        
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);

        /*Call TriaRef function*/
        TriaRef::GetParameterValue(pvalue,&triavalues[0],gauss);

}
/*}}}*/
/*FUNCTION TriaVertexForcing::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/
void TriaVertexForcing::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){

        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);

        /*Call TriaRef function*/
        TriaRef::GetParameterDerivativeValue(p,&triavalues[0],xyz_list,gauss);
}
/*}}}*/
/*FUNCTION TriaVertexForcing::ChangeEnum{{{1*/
void TriaVertexForcing::ChangeEnum(int newenumtype){
        this->enum_type=newenumtype;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::GetParameterAverage{{{1*/
void TriaVertexForcing::GetParameterAverage(double* pvalue){
        
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);

        *pvalue=1./3.*(triavalues[0]+triavalues[1]+triavalues[2]);
}
/*}}}*/

/*Intermediary*/
/*FUNCTION TriaVertexForcing::SquareMin{{{1*/
void TriaVertexForcing::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){

        int i;
        const int numnodes=3;
        double valuescopy[numnodes];
        double squaremin;
        double time;

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&valuescopy[0],time);

        /*Process units if requested: */
        if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type,parameters);

        /*Now, figure out minimum of valuescopy: */
        squaremin=pow(valuescopy[0],2);
        for(i=1;i<numnodes;i++){
                if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
        }
        /*Assign output pointers:*/
        *psquaremin=squaremin;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::InfinityNorm{{{1*/
double TriaVertexForcing::InfinityNorm(void){

        /*Output*/
        double norm=0;
        const int numnodes=3;
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);

        for(int i=0;i<numnodes;i++) if(fabs(triavalues[i])>norm) norm=fabs(triavalues[i]);
        return norm;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::Max{{{1*/
double TriaVertexForcing::Max(void){

        const int numnodes=3;
        double    max;
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);
        
        max=triavalues[0];

        for(int i=1;i<numnodes;i++){
                if(triavalues[i]>max) max=triavalues[i];
        }
        return max;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::MaxAbs{{{1*/
double TriaVertexForcing::MaxAbs(void){

        const int numnodes=3;
        double    max;
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);
        
        
        max=fabs(triavalues[0]);

        for(int i=1;i<numnodes;i++){
                if(fabs(triavalues[i])>max) max=fabs(triavalues[i]);
        }
        return max;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::Min{{{1*/
double TriaVertexForcing::Min(void){

        const int numnodes=3;
        double    min;
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);
        
        min=triavalues[0];

        for(int i=1;i<numnodes;i++){
                if(triavalues[i]<min) min=triavalues[i];
        }
        return min;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::MinAbs{{{1*/
double TriaVertexForcing::MinAbs(void){

        const int numnodes=3;
        double    min;
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);
        
        min=fabs(triavalues[0]);

        for(int i=1;i<numnodes;i++){
                if(fabs(triavalues[i])<min) min=fabs(triavalues[i]);
        }
        return min;
}
/*}}}*/
/*FUNCTION TriaVertexForcing::GetVectorFromInputs{{{1*/
void TriaVertexForcing::GetVectorFromInputs(Vec vector,int* doflist){

        const int numvertices=3;
        double time;
        double triavalues[3];

        /*First, recover current time from parameters: */
        parameters->FindParam(&time,TimeEnum);

        /*Retrieve interpolated values for this time step: */
        this->GetTimeValues(&triavalues[0],time);

        VecSetValues(vector,numvertices,doflist,(const double*)&triavalues[0],INSERT_VALUES);

} /*}}}*/
/*FUNCTION TriaVertexForcing::GetTimeValues{{{1*/
void TriaVertexForcing::GetTimeValues(double* outvalues,double intime){

        int    i,j;
        double deltat;
        double alpha;
        bool   found=false;

        /*Ok, we have the time, go through the timesteps, and figure out which interval we 
         *fall within. Then interpolate the values on this interval: */
        if(intime<this->timesteps[0]){
                /*get values for the first time: */
                outvalues[0]=this->values[3*0+0];
                outvalues[1]=this->values[3*0+1];
                outvalues[2]=this->values[3*0+2];
                found=true;
        }
        else if(intime>this->timesteps[this->numtimesteps-1]){
                /*get values for the last time: */
                outvalues[0]=this->values[3*(this->numtimesteps-1)+0];
                outvalues[1]=this->values[3*(this->numtimesteps-1)+1];
                outvalues[2]=this->values[3*(this->numtimesteps-1)+2];
                found=true;
        }
        else{
                /*Find which interval we fall within: */
                for(i=0;i<this->numtimesteps;i++){
                        if(intime==this->timesteps[i]){
                                /*We are right on one step time: */
                                outvalues[0]=this->values[3*i+0];
                                outvalues[1]=this->values[3*i+1];
                                outvalues[2]=this->values[3*i+2];
                                found=true;
                                break; //we are done with the time interpolation.
                        }
                        else{
                                if(this->timesteps[i]<intime && intime<this->timesteps[i+1]){
                                        /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
                                        deltat=this->timesteps[i+1]-this->timesteps[i];
                                        alpha=(intime-this->timesteps[i])/deltat;
                                        for(j=0;j<3;j++){
                                                outvalues[j]=this->values[3*i+j]+alpha*(this->values[3*(i+1)+j]-this->values[3*i+j]);
                                        }
                                        found=true;
                                        break;
                                }
                                else continue; //keep looking on the next interval
                        }
                }
        }
        if(!found)_error_("did not find time interval on which to interpolate forcing values!");
        
}
/*}}}*/
/*FUNCTION TriaVertexForcing::Configure{{{1*/
void TriaVertexForcing::Configure(Parameters* parameters){
        this->parameters=parameters;
}
/*}}}*/