source: issm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.cpp@ 15300

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

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

#include "../classes.h"
#include "../../shared/shared.h"

/*PentaP1Input constructors and destructor*/
/*FUNCTION PentaP1Input::PentaP1Input(){{{*/
PentaP1Input::PentaP1Input(){
        return;
}
/*}}}*/
/*FUNCTION PentaP1Input::PentaP1Input(int in_enum_type,IssmDouble* values){{{*/
PentaP1Input::PentaP1Input(int in_enum_type,IssmDouble* in_values)
                :PentaRef(1)
{

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

        enum_type=in_enum_type;
        values[0]=in_values[0];
        values[1]=in_values[1];
        values[2]=in_values[2];
        values[3]=in_values[3];
        values[4]=in_values[4];
        values[5]=in_values[5];
}
/*}}}*/
/*FUNCTION PentaP1Input::~PentaP1Input(){{{*/
PentaP1Input::~PentaP1Input(){
        return;
}
/*}}}*/

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

        _printf_("PentaP1Input:\n");
        _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
        _printf_("   values: [" << this->values[0] << " " << this->values[1] << " " << this->values[2] << " " << this->values[3] << " " << this->values[4] << " " << this->values[5] << "]\n");
}
/*}}}*/
/*FUNCTION PentaP1Input::Id{{{*/
int    PentaP1Input::Id(void){ return -1; }
/*}}}*/
/*FUNCTION PentaP1Input::ObjectEnum{{{*/
int PentaP1Input::ObjectEnum(void){

        return PentaP1InputEnum;

}
/*}}}*/

/*PentaP1Input management*/
/*FUNCTION PentaP1Input::copy{{{*/
Object* PentaP1Input::copy() {

        return new PentaP1Input(this->enum_type,this->values);

}
/*}}}*/
/*FUNCTION PentaP1Input::InstanceEnum{{{*/
int PentaP1Input::InstanceEnum(void){

        return this->enum_type;

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

        /*output*/
        TriaInput* outinput=NULL;
        IssmDouble newvalues[3];

        /*Loop over the new indices*/
        for(int i=0;i<3;i++){

                /*Check index value*/
                _assert_(indices[i]>=0 && indices[i]<6);

                /*Assign value to new input*/
                newvalues[i]=this->values[indices[i]];
        }

        /*Create new Tria input*/
        outinput=new TriaInput(this->enum_type,&newvalues[0],P1Enum);

        /*Assign output*/
        return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::SpawnResult{{{*/
ElementResult* PentaP1Input::SpawnResult(int step, IssmDouble time){

        return new PentaP1ElementResult(this->enum_type,this->values,step,time);

}
/*}}}*/

/*Object functions*/
/*FUNCTION PentaP1Input::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/
void PentaP1Input::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){

        /*Call PentaRef function*/
        PentaRef::GetInputValue(pvalue,&values[0],gauss);

}
/*}}}*/
/*FUNCTION PentaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
void PentaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){

        /*Call PentaRef function*/
        PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
}
/*}}}*/
/*FUNCTION PentaP1Input::GetVxStrainRate3d{{{*/
void PentaP1Input::GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
        int i,j;

        const int numnodes=6;
        const int DOFVELOCITY=3;
        IssmDouble B[8][27];
        IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
        IssmDouble velocity[numnodes][DOFVELOCITY];

        /*Get B matrix: */
        GetBStokes(&B[0][0], xyz_list, gauss);
        /*Create a reduced matrix of B to get rid of pressure */
        for (i=0;i<6;i++){
                for (j=0;j<3;j++){
                        B_reduced[i][j]=B[i][j];
                }
                for (j=4;j<7;j++){
                        B_reduced[i][j-1]=B[i][j];
                }
                for (j=8;j<11;j++){
                        B_reduced[i][j-2]=B[i][j];
                }
                for (j=12;j<15;j++){
                        B_reduced[i][j-3]=B[i][j];
                }
                for (j=16;j<19;j++){
                        B_reduced[i][j-4]=B[i][j];
                }
                for (j=20;j<23;j++){
                        B_reduced[i][j-5]=B[i][j];
                }
        }

        /*Here, we are computing the strain rate of (vx,0,0)*/
        for(i=0;i<numnodes;i++){
                velocity[i][0]=this->values[i];
                velocity[i][1]=0.0;
                velocity[i][2]=0.0;
        }
        /*Multiply B by velocity, to get strain rate: */
        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvx,0);

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVyStrainRate3d{{{*/
void PentaP1Input::GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){
        int i,j;

        const int numnodes=6;
        const int DOFVELOCITY=3;
        IssmDouble B[8][27];
        IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
        IssmDouble velocity[numnodes][DOFVELOCITY];

        /*Get B matrix: */
        GetBStokes(&B[0][0], xyz_list, gauss);
        /*Create a reduced matrix of B to get rid of pressure */
        for (i=0;i<6;i++){
                for (j=0;j<3;j++){
                        B_reduced[i][j]=B[i][j];
                }
                for (j=4;j<7;j++){
                        B_reduced[i][j-1]=B[i][j];
                }
                for (j=8;j<11;j++){
                        B_reduced[i][j-2]=B[i][j];
                }
                for (j=12;j<15;j++){
                        B_reduced[i][j-3]=B[i][j];
                }
                for (j=16;j<19;j++){
                        B_reduced[i][j-4]=B[i][j];
                }
                for (j=20;j<23;j++){
                        B_reduced[i][j-5]=B[i][j];
                }
        }

        /*Here, we are computing the strain rate of (0,vy,0)*/
        for(i=0;i<numnodes;i++){
                velocity[i][0]=0.0;
                velocity[i][1]=this->values[i];
                velocity[i][2]=0.0;
        }
        /*Multiply B by velocity, to get strain rate: */
        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvy,0);

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVzStrainRate3d{{{*/
void PentaP1Input::GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss){
        int i,j;

        const int numnodes=6;
        const int DOFVELOCITY=3;
        IssmDouble B[8][27];
        IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
        IssmDouble velocity[numnodes][DOFVELOCITY];

        /*Get B matrix: */
        GetBStokes(&B[0][0], xyz_list, gauss);
        /*Create a reduced matrix of B to get rid of pressure */
        for (i=0;i<6;i++){
                for (j=0;j<3;j++){
                        B_reduced[i][j]=B[i][j];
                }
                for (j=4;j<7;j++){
                        B_reduced[i][j-1]=B[i][j];
                }
                for (j=8;j<11;j++){
                        B_reduced[i][j-2]=B[i][j];
                }
                for (j=12;j<15;j++){
                        B_reduced[i][j-3]=B[i][j];
                }
                for (j=16;j<19;j++){
                        B_reduced[i][j-4]=B[i][j];
                }
                for (j=20;j<23;j++){
                        B_reduced[i][j-5]=B[i][j];
                }
        }

        /*Here, we are computing the strain rate of (0,0,vz)*/
        for(i=0;i<numnodes;i++){
                velocity[i][0]=0.0;
                velocity[i][1]=0.0;
                velocity[i][2]=this->values[i];
        }

        /*Multiply B by velocity, to get strain rate: */
        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvz,0);

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVxStrainRate3dPattyn{{{*/
void PentaP1Input::GetVxStrainRate3dPattyn(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){

        int i;
        const int numnodes=6;
        IssmDouble B[5][NDOF2*numnodes];
        IssmDouble velocity[numnodes][NDOF2];

        /*Get B matrix: */
        GetBPattyn(&B[0][0], xyz_list, gauss);

        /*Here, we are computing the strain rate of (vx,0)*/
        for(i=0;i<numnodes;i++){
                velocity[i][0]=this->values[i];
                velocity[i][1]=0.0;
        }

        /*Multiply B by velocity, to get strain rate: */
        MatrixMultiply( &B[0][0],5,NDOF2*numnodes,0,
                                &velocity[0][0],NDOF2*numnodes,1,0,
                                epsilonvx,0);

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVyStrainRate3dPattyn{{{*/
void PentaP1Input::GetVyStrainRate3dPattyn(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){

        int i;
        const int numnodes=6;
        IssmDouble B[5][NDOF2*numnodes];
        IssmDouble velocity[numnodes][NDOF2];

        /*Get B matrix: */
        GetBPattyn(&B[0][0], xyz_list, gauss);

        /*Here, we are computing the strain rate of (0,vy)*/
        for(i=0;i<numnodes;i++){
                velocity[i][0]=0.0;
                velocity[i][1]=this->values[i];
        }

        /*Multiply B by velocity, to get strain rate: */
        MatrixMultiply( &B[0][0],5,NDOF2*numnodes,0,
                                &velocity[0][0],NDOF2*numnodes,1,0,
                                epsilonvy,0);

}
/*}}}*/
/*FUNCTION PentaP1Input::ChangeEnum{{{*/
void PentaP1Input::ChangeEnum(int newenumtype){
        this->enum_type=newenumtype;
}
/*}}}*/
/*FUNCTION PentaP1Input::GetInputAverage{{{*/
void PentaP1Input::GetInputAverage(IssmDouble* pvalue){
        *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
}
/*}}}*/

/*Intermediary*/
/*FUNCTION PentaP1Input::SquareMin{{{*/
void PentaP1Input::SquareMin(IssmDouble* psquaremin,Parameters* parameters){

        int i;
        const int numnodes=6;
        IssmDouble valuescopy[numnodes];
        IssmDouble squaremin;

        /*First,  copy values, to process units if requested: */
        for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];

        /*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 PentaP1Input::ConstrainMin{{{*/
void PentaP1Input::ConstrainMin(IssmDouble minimum){

        int i;
        const int numnodes=6;

        for(i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
}
/*}}}*/
/*FUNCTION PentaP1Input::InfinityNorm{{{*/
IssmDouble PentaP1Input::InfinityNorm(void){

        /*Output*/
        const int numnodes=6;
        IssmDouble norm=0;

        for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
        return norm;
}
/*}}}*/
/*FUNCTION PentaP1Input::Max{{{*/
IssmDouble PentaP1Input::Max(void){

        const int numnodes=6;
        IssmDouble    max=values[0];

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

        const int numnodes=6;
        IssmDouble    max=fabs(values[0]);

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

        const int numnodes=6;
        IssmDouble    min=values[0];

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

        const int numnodes=6;
        IssmDouble    min=fabs(values[0]);

        for(int i=1;i<numnodes;i++){
                if(fabs(values[i])<min) min=fabs(values[i]);
        }
        return min;
}
/*}}}*/
/*FUNCTION PentaP1Input::Scale{{{*/
void PentaP1Input::Scale(IssmDouble scale_factor){

        int i;
        const int numnodes=6;

        for(i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
}
/*}}}*/
/*FUNCTION PentaP1Input::AXPY{{{*/
void PentaP1Input::AXPY(Input* xinput,IssmDouble scalar){

        int i;
        const int numnodes=6;

        /*xinput is of the same type, so cast it: */

        /*Carry out the AXPY operation depending on type:*/
        switch(xinput->ObjectEnum()){

                case PentaP1InputEnum:{
                        PentaP1Input* cast_input=(PentaP1Input*)xinput;
                        for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
                        return;
                case ControlInputEnum:{
                        ControlInput* cont_input=(ControlInput*)xinput;
                        if(cont_input->values->ObjectEnum()!=PentaP1InputEnum) _error_("not supported yet");
                        PentaP1Input* cast_input=(PentaP1Input*)cont_input->values;
                        for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
                        return;
                default:
                        _error_("not implemented yet");
        }

}
/*}}}*/
/*FUNCTION PentaP1Input::Constrain{{{*/
void PentaP1Input::Constrain(IssmDouble cm_min, IssmDouble cm_max){

        int i;
        const int numnodes=6;

        if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
        if(!xIsNan<IssmDouble>(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;

}
/*}}}*/
/*FUNCTION PentaP1Input::Extrude{{{*/
void PentaP1Input::Extrude(void){

        int i;

        /*First 3 values copied on 3 last values*/
        for(i=0;i<3;i++) this->values[3+i]=this->values[i];
}
/*}}}*/
/*FUNCTION PentaP1Input::VerticallyIntegrate{{{*/
void PentaP1Input::VerticallyIntegrate(Input* thickness_input){

        IssmDouble thickness;

        /*Check that input provided is a thickness*/
        if (thickness_input->InstanceEnum()!=ThicknessEnum) _error_("Input provided is not a Thickness (enum_type is " << EnumToStringx(thickness_input->InstanceEnum()) << ")");

        /*vertically integrate depending on type:*/
        switch(thickness_input->ObjectEnum()){

                case PentaP1InputEnum:{
                        GaussPenta *gauss=new GaussPenta();
                        for(int iv=0;iv<3;iv++){
                                gauss->GaussVertex(iv);
                                thickness_input->GetInputValue(&thickness,gauss);
                                this->values[iv]=0.5*(this->values[iv]+this->values[iv+3]) * thickness;
                                this->values[iv+3]=this->values[iv];
                        }
                        delete gauss;
                        return; }

                default:
                        _error_("not implemented yet");
        }
}
/*}}}*/
/*FUNCTION PentaP1Input::PointwiseDivide{{{*/
Input* PentaP1Input::PointwiseDivide(Input* inputB){

        /*Ouput*/
        PentaP1Input* outinput=NULL;

        /*Intermediaries*/
        PentaP1Input *xinputB  = NULL;
        const int     numnodes = 6;
        IssmDouble    AdotBvalues[numnodes];

        /*Check that inputB is of the same type*/
        if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
        xinputB=(PentaP1Input*)inputB;

        /*Create point wise sum*/
        for(int i=0;i<numnodes;i++){
                _assert_(xinputB->values[i]!=0);
                AdotBvalues[i]=this->values[i]/xinputB->values[i];
        }

        /*Create new Penta vertex input (copy of current input)*/
        outinput=new PentaP1Input(this->enum_type,&AdotBvalues[0]);

        /*Return output pointer*/
        return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::PointwiseMin{{{*/
Input* PentaP1Input::PointwiseMin(Input* inputB){

        /*Ouput*/
        PentaP1Input* outinput=NULL;

        /*Intermediaries*/
        int               i;
        PentaP1Input *xinputB     = NULL;
        const int         numnodes    = 6;
        IssmDouble            minvalues[numnodes];

        /*Check that inputB is of the same type*/
        if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
        xinputB=(PentaP1Input*)inputB;

        /*Create point wise min*/
        for(i=0;i<numnodes;i++){
                if(this->values[i] > xinputB->values[i]) minvalues[i]=xinputB->values[i];
                else minvalues[i]=this->values[i];
        }

        /*Create new Penta vertex input (copy of current input)*/
        outinput=new PentaP1Input(this->enum_type,&minvalues[0]);

        /*Return output pointer*/
        return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::PointwiseMax{{{*/
Input* PentaP1Input::PointwiseMax(Input* inputB){

        /*Ouput*/
        PentaP1Input* outinput=NULL;

        /*Intermediaries*/
        int               i;
        PentaP1Input *xinputB     = NULL;
        const int         numnodes    = 6;
        IssmDouble            maxvalues[numnodes];

        /*Check that inputB is of the same type*/
        if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
        xinputB=(PentaP1Input*)inputB;

        /*Create point wise max*/
        for(i=0;i<numnodes;i++){
                if(this->values[i] < xinputB->values[i]) maxvalues[i]=xinputB->values[i];
                else maxvalues[i]=this->values[i];
        }

        /*Create new Penta vertex input (copy of current input)*/
        outinput=new PentaP1Input(this->enum_type,&maxvalues[0]);

        /*Return output pointer*/
        return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVectorFromInputs{{{*/
void PentaP1Input::GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist){

        const int numvertices=6;
        vector->SetValues(numvertices,doflist,this->values,INS_VAL);

} /*}}}*/
/*FUNCTION PentaP1Input::Configure{{{*/
void PentaP1Input::Configure(Parameters* parameters){
        /*do nothing: */
}
/*}}}*/