source: issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp@ 4899

/*!\file PentaVertexInput.c
 * \brief: implementation of the PentaVertexInput 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"

/*PentaVertexInput constructors and destructor*/
/*FUNCTION PentaVertexInput::PentaVertexInput(){{{1*/
PentaVertexInput::PentaVertexInput(){
        return;
}
/*}}}*/
/*FUNCTION PentaVertexInput::PentaVertexInput(int in_enum_type,double* values){{{1*/
PentaVertexInput::PentaVertexInput(int in_enum_type,double* 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 PentaVertexInput::~PentaVertexInput(){{{1*/
PentaVertexInput::~PentaVertexInput(){
        return;
}
/*}}}*/

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

        printf("PentaVertexInput:\n");
        printf("   enum: %i (%s)\n",this->enum_type,EnumAsString(this->enum_type));
        printf("   values: [%g %g %g %g %g %g]\n",this->values[0],this->values[1],this->values[2],this->values[3],this->values[4],this->values[5]);
}
/*}}}*/
/*FUNCTION PentaVertexInput::Id{{{1*/
int    PentaVertexInput::Id(void){ return -1; }
/*}}}*/
/*FUNCTION PentaVertexInput::MyRank{{{1*/
int    PentaVertexInput::MyRank(void){ 
        extern int my_rank;
        return my_rank; 
}
/*}}}*/
/*FUNCTION PentaVertexInput::Marshall{{{1*/
void  PentaVertexInput::Marshall(char** pmarshalled_dataset){

        char* marshalled_dataset=NULL;
        int   enum_value=0;

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

        /*get enum value of PentaVertexInput: */
        enum_value=PentaVertexInputEnum;
        
        /*marshall enum: */
        memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
        
        /*marshall PentaVertexInput data: */
        memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
        memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values);

        *pmarshalled_dataset=marshalled_dataset;
}
/*}}}*/
/*FUNCTION PentaVertexInput::MarshallSize{{{1*/
int   PentaVertexInput::MarshallSize(){
        
        return sizeof(values)+
                +sizeof(enum_type)+
                +sizeof(int); //sizeof(int) for enum value
}
/*}}}*/
/*FUNCTION PentaVertexInput::Demarshall{{{1*/
void  PentaVertexInput::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(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values);

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

        return PentaVertexInputEnum;

}
/*}}}*/
        
/*PentaVertexInput management*/
/*FUNCTION PentaVertexInput::copy{{{1*/
Object* PentaVertexInput::copy() {
        
        return new PentaVertexInput(this->enum_type,this->values);

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

        return this->enum_type;

}
/*}}}*/
/*FUNCTION PentaVertexInput::SpawnBeamInput{{{1*/
Input* PentaVertexInput::SpawnBeamInput(int* indices){

        /*output*/
        BeamVertexInput* outinput=NULL;
        double newvalues[2];

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

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

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

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

        /*Assign output*/
        return outinput;

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

        /*output*/
        TriaVertexInput* outinput=NULL;
        double newvalues[3];

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

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

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

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

        /*Assign output*/
        return outinput;

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

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

}
/*}}}*/

/*Object functions*/
/*FUNCTION PentaVertexInput::GetParameterValue(bool* pvalue) {{{1*/
void PentaVertexInput::GetParameterValue(bool* pvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(int* pvalue){{{1*/
void PentaVertexInput::GetParameterValue(int* pvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
        /*P1 interpolation on Gauss point*/

        /*intermediary*/
        double l1l6[6];

        /*nodal functions: */
        GetNodalFunctionsP1(&l1l6[0],gauss);

        /*Assign output pointers:*/
        *pvalue=l1l6[0]*values[0]+l1l6[1]*values[1]+l1l6[2]*values[2]+l1l6[3]*values[3]+l1l6[4]*values[4]+l1l6[5]*values[5];

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValues{{{1*/
void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){
        /*It is assumed that output values has been correctly allocated*/

        int i,j;
        double gauss[4];

        for (i=0;i<numgauss;i++){

                /*Get current Gauss point coordinates*/
                for (j=0;j<4;j++) gauss[j]=gauss_pointers[i*4+j];

                /*Assign parameter value*/
                GetParameterValue(&values[i],&gauss[0]);
        }
}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterDerivativeValue{{{1*/
void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
        /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
         *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
         *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
         *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
         *
         *   p is a vector of size 3x1 already allocated.
         */

        const int NDOF3=3;
        const int numgrids=6;
        double dh1dh6[NDOF3][numgrids];

        /*Get nodal funnctions derivatives in actual coordinate system: */
        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);

        p[0]=this->values[0]*dh1dh6[0][0]+this->values[1]*dh1dh6[0][1]+this->values[2]*dh1dh6[0][2]+this->values[3]*dh1dh6[0][3]+this->values[4]*dh1dh6[0][4]+this->values[5]*dh1dh6[0][5];
        p[1]=this->values[0]*dh1dh6[1][0]+this->values[1]*dh1dh6[1][1]+this->values[2]*dh1dh6[1][2]+this->values[3]*dh1dh6[1][3]+this->values[4]*dh1dh6[1][4]+this->values[5]*dh1dh6[1][5];
        p[2]=this->values[0]*dh1dh6[2][0]+this->values[1]*dh1dh6[2][1]+this->values[2]*dh1dh6[2][2]+this->values[3]*dh1dh6[2][3]+this->values[4]*dh1dh6[2][4]+this->values[5]*dh1dh6[2][5];

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
        int i,j;

        const int numgrids=6;
        const int DOFVELOCITY=3;
        double B[8][27];
        double B_reduced[6][DOFVELOCITY*numgrids];
        double velocity[numgrids][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<numgrids;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*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
        int i,j;

        const int numgrids=6;
        const int DOFVELOCITY=3;
        double B[8][27];
        double B_reduced[6][DOFVELOCITY*numgrids];
        double velocity[numgrids][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<numgrids;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*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
        int i,j;

        const int numgrids=6;
        const int DOFVELOCITY=3;
        double B[8][27];
        double B_reduced[6][DOFVELOCITY*numgrids];
        double velocity[numgrids][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<numgrids;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*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){

        int i;
        const int numgrids=6;
        const int NDOF2=2;
        double B[5][NDOF2*numgrids];
        double velocity[numgrids][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<numgrids;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*numgrids,0,
                                &velocity[0][0],NDOF2*numgrids,1,0,
                                epsilonvx,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){

        int i;
        const int numgrids=6;
        const int NDOF2=2;
        double B[5][NDOF2*numgrids];
        double velocity[numgrids][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<numgrids;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*numgrids,0,
                                &velocity[0][0],NDOF2*numgrids,1,0,
                                epsilonvy,0);

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

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

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

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

        /*Process units if requested: */
        if(process_units)NodalValuesUnitConversion(&valuescopy[0],numnodes,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 PentaVertexInput::Scale{{{1*/
void PentaVertexInput::Scale(double scale_factor){
        
        int i;
        const int numgrids=6;

        for(i=0;i<numgrids;i++)values[i]=values[i]*scale_factor;
}
/*}}}*/
/*FUNCTION PentaVertexInput::AXPY{{{1*/
void PentaVertexInput::AXPY(Input* xinput,double scalar){

        int i;
        const int numgrids=6;
        PentaVertexInput*  xpentavertexinput=NULL;

        /*xinput is of the same type, so cast it: */
        xpentavertexinput=(PentaVertexInput*)xinput;

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

                case PentaVertexInputEnum:
                        for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
                        return;

                default:
                        ISSMERROR("not implemented yet");
        }

}
/*}}}*/
/*FUNCTION PentaVertexInput::Constrain{{{1*/
void PentaVertexInput::Constrain(double cm_min, double cm_max){

        int i;
        const int numgrids=6;
                
        if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
        if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;

}
/*}}}*/
/*FUNCTION PentaVertexInput::Extrude{{{1*/
void PentaVertexInput::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 PentaVertexInput::VerticallyIntegrate{{{1*/
void PentaVertexInput::VerticallyIntegrate(Input* thickness_input){

        /*Intermediaries*/
        int i;
        const int  numgrids = 6;
        int        num_thickness_values;
        double    *thickness_values = NULL;

        /*Check that input provided is a thickness*/
        if (thickness_input->EnumType()!=ThicknessEnum) ISSMERROR("Input provided is not a Thickness (enum_type is %s)",EnumAsString(thickness_input->EnumType()));

        /*Get Thickness value pointer*/
        thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);

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

                case PentaVertexInputEnum:
                        for(i=0;i<3;i++){
                                this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
                                this->values[i+3]=this->values[i];
                        }
                        return;

                default:
                        ISSMERROR("not implemented yet");
        }
}
/*}}}*/
/*FUNCTION PentaVertexInput::PointwiseDivide{{{1*/
Input* PentaVertexInput::PointwiseDivide(Input* inputB){

        /*Ouput*/
        PentaVertexInput* outinput=NULL;

        /*Intermediaries*/
        int               i;
        PentaVertexInput *xinputB     = NULL;
        int               B_numvalues;
        double           *B_values    = NULL;
        const int         numgrids    = 6;
        double            AdotBvalues[numgrids];

        /*Check that inputB is of the same type*/
        if (inputB->Enum()!=PentaVertexInputEnum) ISSMERROR("Operation not permitted because inputB is of type %s",EnumAsString(inputB->Enum()));
        xinputB=(PentaVertexInput*)inputB;

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

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

        /*Return output pointer*/
        return outinput;

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVectorFromInputs{{{1*/
void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){

        const int numvertices=6;
        VecSetValues(vector,numvertices,doflist,(const double*)this->values,INSERT_VALUES);

} /*}}}*/
/*FUNCTION PentaVertexInput::GetValuesPtr{{{1*/
void PentaVertexInput::GetValuesPtr(double** pvalues,int* pnum_values){

        *pvalues=this->values;
        *pnum_values=6;

}
/*}}}*/