/*!\file PentaP1Input.c * \brief: implementation of the PentaP1Input object */ #ifdef HAVE_CONFIG_H #include #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif #include #include #include "../objects.h" #include "../../EnumDefinitions/EnumDefinitions.h" #include "../../shared/shared.h" #include "../../Container/Container.h" #include "../../include/include.h" /*PentaP1Input constructors and destructor*/ /*FUNCTION PentaP1Input::PentaP1Input(){{{*/ PentaP1Input::PentaP1Input(){ return; } /*}}}*/ /*FUNCTION PentaP1Input::PentaP1Input(int in_enum_type,IssmPDouble* values){{{*/ PentaP1Input::PentaP1Input(int in_enum_type,IssmPDouble* 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){ _printLine_("PentaP1Input:"); _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); _printLine_(" values: [" << this->values[0] << " " << this->values[1] << " " << this->values[2] << " " << this->values[3] << " " << this->values[4] << " " << this->values[5] << "]"); } /*}}}*/ /*FUNCTION PentaP1Input::Id{{{*/ int PentaP1Input::Id(void){ return -1; } /*}}}*/ /*FUNCTION PentaP1Input::MyRank{{{*/ int PentaP1Input::MyRank(void){ extern int my_rank; return my_rank; } /*}}}*/ /*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*/ TriaP1Input* outinput=NULL; IssmPDouble 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 TriaP1Input(this->enum_type,&newvalues[0]); /*Assign output*/ return outinput; } /*}}}*/ /*FUNCTION PentaP1Input::SpawnResult{{{*/ ElementResult* PentaP1Input::SpawnResult(int step, IssmPDouble time){ return new PentaP1ElementResult(this->enum_type,this->values,step,time); } /*}}}*/ /*Object functions*/ /*FUNCTION PentaP1Input::GetInputValue(IssmPDouble* pvalue,GaussPenta* gauss){{{*/ void PentaP1Input::GetInputValue(IssmPDouble* pvalue,GaussPenta* gauss){ /*Call PentaRef function*/ PentaRef::GetInputValue(pvalue,&values[0],gauss); } /*}}}*/ /*FUNCTION PentaP1Input::GetInputDerivativeValue(IssmPDouble* p, IssmPDouble* xyz_list, GaussPenta* gauss){{{*/ void PentaP1Input::GetInputDerivativeValue(IssmPDouble* p, IssmPDouble* xyz_list, GaussPenta* gauss){ /*Call PentaRef function*/ PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss); } /*}}}*/ /*FUNCTION PentaP1Input::GetVxStrainRate3d{{{*/ void PentaP1Input::GetVxStrainRate3d(IssmPDouble* epsilonvx,IssmPDouble* xyz_list, GaussPenta* gauss){ int i,j; const int numnodes=6; const int DOFVELOCITY=3; IssmPDouble B[8][27]; IssmPDouble B_reduced[6][DOFVELOCITY*numnodes]; IssmPDouble 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;ivalues[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(IssmPDouble* epsilonvy,IssmPDouble* xyz_list, GaussPenta* gauss){ int i,j; const int numnodes=6; const int DOFVELOCITY=3; IssmPDouble B[8][27]; IssmPDouble B_reduced[6][DOFVELOCITY*numnodes]; IssmPDouble 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;ivalues[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(IssmPDouble* epsilonvz,IssmPDouble* xyz_list, GaussPenta* gauss){ int i,j; const int numnodes=6; const int DOFVELOCITY=3; IssmPDouble B[8][27]; IssmPDouble B_reduced[6][DOFVELOCITY*numnodes]; IssmPDouble 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;ivalues[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(IssmPDouble* epsilonvx,IssmPDouble* xyz_list, GaussPenta* gauss){ int i; const int numnodes=6; IssmPDouble B[5][NDOF2*numnodes]; IssmPDouble 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;ivalues[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(IssmPDouble* epsilonvy,IssmPDouble* xyz_list, GaussPenta* gauss){ int i; const int numnodes=6; IssmPDouble B[5][NDOF2*numnodes]; IssmPDouble 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;ivalues[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(IssmPDouble* pvalue){ *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]); } /*}}}*/ /*Intermediary*/ /*FUNCTION PentaP1Input::SquareMin{{{*/ void PentaP1Input::SquareMin(IssmPDouble* psquaremin, bool process_units,Parameters* parameters){ int i; const int numnodes=6; IssmPDouble valuescopy[numnodes]; IssmPDouble squaremin; /*First, copy values, to process units if requested: */ for(i=0;ivalues[i]; /*Process units if requested: */ if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type); /*Now, figure out minimum of valuescopy: */ squaremin=pow(valuescopy[0],2); for(i=1;inorm) norm=fabs(values[i]); return norm; } /*}}}*/ /*FUNCTION PentaP1Input::Max{{{*/ IssmPDouble PentaP1Input::Max(void){ const int numnodes=6; IssmPDouble max=values[0]; for(int i=1;imax) max=values[i]; } return max; } /*}}}*/ /*FUNCTION PentaP1Input::MaxAbs{{{*/ IssmPDouble PentaP1Input::MaxAbs(void){ const int numnodes=6; IssmPDouble max=fabs(values[0]); for(int i=1;imax) max=fabs(values[i]); } return max; } /*}}}*/ /*FUNCTION PentaP1Input::Min{{{*/ IssmPDouble PentaP1Input::Min(void){ const int numnodes=6; IssmPDouble min=values[0]; for(int i=1;iObjectEnum()){ case PentaP1InputEnum:{ PentaP1Input* cast_input=(PentaP1Input*)xinput; for(i=0;ivalues[i]=this->values[i]+scalar*(cast_input->values[i]);} return; case ControlInputEnum:{ ControlInput* cont_input=(ControlInput*)xinput; if(cont_input->values->ObjectEnum()!=PentaP1InputEnum) _error2_("not supported yet"); PentaP1Input* cast_input=(PentaP1Input*)cont_input->values; for(i=0;ivalues[i]=this->values[i]+scalar*(cast_input->values[i]);} return; default: _error2_("not implemented yet"); } } /*}}}*/ /*FUNCTION PentaP1Input::Constrain{{{*/ void PentaP1Input::Constrain(IssmPDouble cm_min, IssmPDouble cm_max){ int i; const int numnodes=6; if(!xIsNan(cm_min)) for(i=0;ivalues[i]values[i]=cm_min; if(!xIsNan(cm_max)) for(i=0;ivalues[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){ /*Intermediaries*/ int i; const int numnodes = 6; int num_thickness_values; IssmPDouble *thickness_values = NULL; /*Check that input provided is a thickness*/ if (thickness_input->InstanceEnum()!=ThicknessEnum) _error2_("Input provided is not a Thickness (enum_type is " << EnumToStringx(thickness_input->InstanceEnum()) << ")"); /*Get Thickness value pointer*/ thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values); /*vertically integrate depending on type:*/ switch(thickness_input->ObjectEnum()){ case PentaP1InputEnum: 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: _error2_("not implemented yet"); } } /*}}}*/ /*FUNCTION PentaP1Input::PointwiseDivide{{{*/ Input* PentaP1Input::PointwiseDivide(Input* inputB){ /*Ouput*/ PentaP1Input* outinput=NULL; /*Intermediaries*/ int i; PentaP1Input *xinputB = NULL; int B_numvalues; const int numnodes = 6; IssmPDouble AdotBvalues[numnodes]; /*Check that inputB is of the same type*/ if (inputB->ObjectEnum()!=PentaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); xinputB=(PentaP1Input*)inputB; /*Create point wise sum*/ for(i=0;ivalues[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; int B_numvalues; const int numnodes = 6; IssmPDouble minvalues[numnodes]; /*Check that inputB is of the same type*/ if (inputB->ObjectEnum()!=PentaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); xinputB=(PentaP1Input*)inputB; /*Create point wise min*/ for(i=0;ivalues[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; int B_numvalues; const int numnodes = 6; IssmPDouble maxvalues[numnodes]; /*Check that inputB is of the same type*/ if (inputB->ObjectEnum()!=PentaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); xinputB=(PentaP1Input*)inputB; /*Create point wise max*/ for(i=0;ivalues[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* vector,int* doflist){ const int numvertices=6; vector->SetValues(numvertices,doflist,this->values,INS_VAL); } /*}}}*/ /*FUNCTION PentaP1Input::GetValuesPtr{{{*/ void PentaP1Input::GetValuesPtr(IssmPDouble** pvalues,int* pnum_values){ *pvalues=this->values; *pnum_values=6; } /*}}}*/ /*FUNCTION PentaP1Input::Configure{{{*/ void PentaP1Input::Configure(Parameters* parameters){ /*do nothing: */ } /*}}}*/