Index: /issm/trunk/src/c/modules/InputDepthAveragex/InputDepthAveragex.cpp =================================================================== --- /issm/trunk/src/c/modules/InputDepthAveragex/InputDepthAveragex.cpp (revision 4470) +++ /issm/trunk/src/c/modules/InputDepthAveragex/InputDepthAveragex.cpp (revision 4471) @@ -25,6 +25,5 @@ } - /*Then extrude vertically the new inputs*/ - InputExtrudex( elements,nodes,vertices,loads,materials,parameters,enum_type); + InputExtrudex( elements,nodes,vertices,loads,materials,parameters,average_enum_type); } Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4470) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4471) @@ -1080,5 +1080,4 @@ /*FUNCTION Penta::InputDepthAverageAtBase{{{1*/ void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){ - ISSMERROR("Not implemented yet (see Node::FieldDepthAverageAtBase)"); /*Intermediaries*/ @@ -1090,5 +1089,6 @@ Penta* penta=NULL; Input* original_input=NULL; - Input* integrated_input=NULL; + Input* element_integrated_input=NULL; + Input* total_integrated_input=NULL; Input* element_thickness_input=NULL; Input* total_thickness_input=NULL; @@ -1102,50 +1102,60 @@ inputs->GetParameterValue(&onbed,ElementOnBedEnum); - /*Are we on the base, not on the surface?:*/ - if(onbed==1){ - - /*OK, we are on bed. Initialize inputs*/ - total_thickness_input=new PentaVertexInput(ThicknessEnum,zeros_list); - depth_averaged_input =new PentaVertexInput(average_enum_type,zeros_list); - - /*Now follow all the upper element from the base to the surface to integrate the input*/ - penta=this; - for(;;){ - - /*Step1: Get original input (to be depth avegaged): */ - original_input=(Input*)penta->inputs->GetInput(enum_type); - if(!original_input) ISSMERROR("%s%s"," could not find input with enum:",EnumAsString(enum_type)); - - /*Step2: Create element thickness input*/ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); - for(i=0;i<3;i++){ - Helem_list[i]=xyz_list[i+3][2]-xyz_list[i][2]; - Helem_list[i+3]=Helem_list[i]; - } - element_thickness_input=new PentaVertexInput(ThicknessEnum,Helem_list); - - /*Step3: Vertically integrate original input and update totalthickness_input*/ - //integrated_input=original_input->VerticallyIntegrate(element_thickness_input); //TO BE ADDED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - total_thickness_input->AXPY(element_thickness_input,1.0); - - /*Stop if we have reached the surface, else, take upper penta*/ - if (penta->IsOnSurface()){ - break; - } - else{ - /* get upper Penta*/ - penta=penta->GetUpperElement(); - ISSMASSERT(penta->Id()!=this->id); - } + /*Are we on the base? If not, return*/ + if(!onbed) return; + + /*OK, we are on bed. Initialize global inputs as 0*/ + total_thickness_input =new PentaVertexInput(ThicknessEnum,zeros_list); + total_integrated_input=new PentaVertexInput(average_enum_type,zeros_list); + + /*Now follow all the upper element from the base to the surface to integrate the input*/ + penta=this; + for(;;){ + + /*Step1: Get original input (to be depth avegaged): */ + original_input=(Input*)penta->inputs->GetInput(enum_type); + if(!original_input) ISSMERROR("could not find input with enum %s",EnumAsString(enum_type)); + + /*Step2: Create element thickness input*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); + for(i=0;i<3;i++){ + Helem_list[i]=xyz_list[i+3][2]-xyz_list[i][2]; + Helem_list[i+3]=Helem_list[i]; } - - /*OK, now we only need to divide the depth integrated input by the total thickness!*/ - //depth_averaged_input=integrated_input->InputPointwiseDivide(total_thickness_input); //TO BE ADDED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - depth_averaged_input->ChangeEnum(average_enum_type); - - /*Finally, add to inputs*/ - this->inputs->AddInput((Input*)depth_averaged_input); - } - + element_thickness_input=new PentaVertexInput(ThicknessEnum,Helem_list); + + /*Step3: Vertically integrate A COPY of the original*/ + element_integrated_input=(Input*)original_input->copy(); + element_integrated_input->VerticallyIntegrate(element_thickness_input); + + /*Add contributions to global inputs*/ + total_integrated_input->AXPY(element_integrated_input,1.0); + total_thickness_input ->AXPY(element_thickness_input,1.0); + + /*Clean up*/ + delete element_thickness_input; + delete element_integrated_input; + + /*Stop if we have reached the surface, else, take upper penta*/ + if (penta->IsOnSurface()){ + break; + } + else{ + /* get upper Penta*/ + penta=penta->GetUpperElement(); + ISSMASSERT(penta->Id()!=this->id); + } + } + + /*OK, now we only need to divide the depth integrated input by the total thickness!*/ + depth_averaged_input=total_integrated_input->PointwiseDivide(total_thickness_input); + depth_averaged_input->ChangeEnum(average_enum_type); + + /*Clean up*/ + delete total_thickness_input; + delete total_integrated_input; + + /*Finally, add to inputs*/ + this->inputs->AddInput((Input*)depth_averaged_input); return; } Index: /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h (revision 4471) @@ -41,4 +41,5 @@ Input* SpawnBeamInput(int* indices); Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");}; ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -70,4 +71,5 @@ void Constrain(double cm_min, double cm_max); void Extrude(void){ISSMERROR("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");}; void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values); Index: /issm/trunk/src/c/objects/Inputs/BoolInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 4471) @@ -41,4 +41,5 @@ Input* SpawnBeamInput(int* indices); Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");}; ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -70,4 +71,5 @@ void Constrain(double cm_min, double cm_max); void Extrude(void){ISSMERROR("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");}; void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values); Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 4471) @@ -40,4 +40,5 @@ Input* SpawnBeamInput(int* indices); Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");}; ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -69,4 +70,5 @@ void Constrain(double cm_min, double cm_max); void Extrude(void){ISSMERROR("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");}; void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values); Index: /issm/trunk/src/c/objects/Inputs/Input.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/Input.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/Input.h (revision 4471) @@ -47,4 +47,5 @@ virtual Input* SpawnBeamInput(int* indices)=0; virtual Input* SpawnTriaInput(int* indices)=0; + virtual Input* PointwiseDivide(Input* inputB)=0; virtual ElementResult* SpawnResult(int step, double time)=0; virtual void SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0; @@ -52,4 +53,5 @@ virtual void AXPY(Input* xinput,double scalar)=0; virtual void Constrain(double cm_min, double cm_max)=0; + virtual void VerticallyIntegrate(Input* thickness_input)=0; virtual void Extrude()=0; virtual void GetVectorFromInputs(Vec vector,int* doflist)=0; Index: /issm/trunk/src/c/objects/Inputs/IntInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 4471) @@ -41,4 +41,5 @@ Input* SpawnBeamInput(int* indices); Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");}; ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -70,4 +71,5 @@ void Constrain(double cm_min, double cm_max); void Extrude(void){ISSMERROR("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");}; void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values); Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 4471) @@ -878,5 +878,5 @@ } /*}}}*/ -/*FUNCTION PentaVertexInput::SquareMin(double* psquaremin, bool process_units){{{1*/ +/*FUNCTION PentaVertexInput::SquareMin{{{1*/ void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){ @@ -901,5 +901,5 @@ } /*}}}*/ -/*FUNCTION PentaVertexInput::Scale(double scale_factor){{{1*/ +/*FUNCTION PentaVertexInput::Scale{{{1*/ void PentaVertexInput::Scale(double scale_factor){ @@ -910,5 +910,5 @@ } /*}}}*/ -/*FUNCTION PentaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/ +/*FUNCTION PentaVertexInput::AXPY{{{1*/ void PentaVertexInput::AXPY(Input* xinput,double scalar){ @@ -933,5 +933,5 @@ } /*}}}*/ -/*FUNCTION PentaVertexInput::Constrain(double cm_min, double cm_max){{{1*/ +/*FUNCTION PentaVertexInput::Constrain{{{1*/ void PentaVertexInput::Constrain(double cm_min, double cm_max){ @@ -944,5 +944,5 @@ } /*}}}*/ -/*FUNCTION PentaVertexInput::Extrude(double cm_min, double cm_max){{{1*/ +/*FUNCTION PentaVertexInput::Extrude{{{1*/ void PentaVertexInput::Extrude(void){ @@ -954,4 +954,66 @@ } /*}}}*/ +/*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;ivalues[i]!=0); + AdotBvalues[i]=this->values[i]/xinputB->values[i]; + } + + /*Create new Sing input (copy of current input)*/ + outinput=new PentaVertexInput(this->enum_type,&AdotBvalues[0]); + + /*Return output pointer*/ + return outinput; + +} +/*}}}*/ /*FUNCTION PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/ void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){ Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 4471) @@ -40,4 +40,5 @@ Input* SpawnBeamInput(int* indices); Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB); ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -79,4 +80,5 @@ void Constrain(double cm_min, double cm_max); void Extrude(void); + void VerticallyIntegrate(Input* thickness_input); void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values); Index: /issm/trunk/src/c/objects/Inputs/SingVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/SingVertexInput.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/SingVertexInput.h (revision 4471) @@ -40,4 +40,5 @@ Input* SpawnBeamInput(int* indices); Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");}; ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -69,4 +70,5 @@ void Constrain(double cm_min, double cm_max); void Extrude(void){ISSMERROR("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");}; void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values); Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 4470) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 4471) @@ -41,4 +41,5 @@ Input* SpawnBeamInput(int* indices); Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");}; ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -77,4 +78,5 @@ void Constrain(double cm_min, double cm_max); void Extrude(void){ISSMERROR("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");}; void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values);