Index: /issm/trunk/src/c/modules/InputDepthAveragex/InputDepthAveragex.cpp =================================================================== --- /issm/trunk/src/c/modules/InputDepthAveragex/InputDepthAveragex.cpp (revision 4930) +++ /issm/trunk/src/c/modules/InputDepthAveragex/InputDepthAveragex.cpp (revision 4931) @@ -19,5 +19,5 @@ for (i=0;iSize();i++){ element=(Element*)elements->GetObjectByOffset(i); - element->InputDepthAverageAtBase(enum_type,average_enum_type); + element->InputDepthAverageAtBase(enum_type,average_enum_type,ElementsEnum); } Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 4930) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 4931) @@ -43,5 +43,5 @@ virtual double CostFunction(void)=0; virtual double SurfaceArea(void)=0; - virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type)=0; + virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0; virtual void ComputeBasalStress(Vec sigma_b)=0; virtual void ComputePressure(Vec p_g)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4930) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4931) @@ -14,5 +14,4 @@ #include #include "../objects.h" -#include "../../EnumDefinitions/EnumDefinitions.h" #include "../../shared/shared.h" #include "../../include/include.h" @@ -1028,10 +1027,10 @@ /*}}}*/ /*FUNCTION Penta::InputDepthAverageAtBase{{{1*/ -void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){ +void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ /*Intermediaries*/ const int numvertices=6; bool onbed; - int i; + int step,i; Penta* penta=NULL; @@ -1055,13 +1054,28 @@ /*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; + step =0; for(;;){ /*Step1: Get original input (to be depth avegaged): */ - original_input=(Input*)penta->inputs->GetInput(enum_type); + if (object_enum==ElementsEnum) + original_input=(Input*)penta->inputs->GetInput(enum_type); + else if (object_enum==MaterialsEnum) + original_input=(Input*)penta->matice->inputs->GetInput(enum_type); + else + ISSMERROR("object %s not supported yet",EnumAsString(object_enum)); if(!original_input) ISSMERROR("could not find input with enum %s",EnumAsString(enum_type)); + + /*If first time, initialize total_integrated_input*/ + if (step==0){ + if (original_input->Enum()==PentaVertexInputEnum) + total_integrated_input=new PentaVertexInput(average_enum_type,zeros_list); + else if (original_input->Enum()==DoubleInputEnum) + total_integrated_input=new DoubleInput(average_enum_type,0.0); + else + ISSMERROR("object %s not supported yet",EnumAsString(original_input->Enum())); + } /*Step2: Create element thickness input*/ @@ -1086,12 +1100,12 @@ /*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); - } + if (penta->IsOnSurface()) break; + + /* get upper Penta*/ + penta=penta->GetUpperElement(); + ISSMASSERT(penta->Id()!=this->id); + + /*increase couter*/ + step++; } @@ -1105,6 +1119,10 @@ /*Finally, add to inputs*/ - this->inputs->AddInput((Input*)depth_averaged_input); - return; + if (object_enum==ElementsEnum) + this->inputs->AddInput((Input*)depth_averaged_input); + else if (object_enum==MaterialsEnum) + this->matice->inputs->AddInput((Input*)depth_averaged_input); + else + ISSMERROR("object %s not supported yet",EnumAsString(object_enum)); } /*}}}*/ @@ -2130,7 +2148,16 @@ /*This element should be collapsed into a tria element at its base. Create this tria element, *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ + + /*Depth Averaging B*/ + this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum); + + /*Call Tria function*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. tria->CreateKMatrix(Kgg); delete tria; + + /*Delete B averaged*/ + this->matice->inputs->DeleteInput(RheologyB2dEnum); + return; } @@ -3941,6 +3968,6 @@ /*Delete Vx and Vy averaged*/ - this-inputs->DeleteInput(VxAverageEnum); - this-inputs->DeleteInput(VyAverageEnum); + this->inputs->DeleteInput(VxAverageEnum); + this->inputs->DeleteInput(VyAverageEnum); return; @@ -5228,7 +5255,7 @@ this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,values)); - /*Also update the rheology WARNING: ONLY FOR PATTYN and STOKES FOR NOW*/ + /*Also update the rheology*/ inputs->GetParameterValue(&collapse,CollapseEnum); - if (!collapse) this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); + this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); } Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 4930) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 4931) @@ -22,4 +22,5 @@ #include "../../include/include.h" #include "../../shared/Exceptions/exceptions.h" +#include "../../EnumDefinitions/EnumDefinitions.h" /*}}}*/ @@ -88,5 +89,5 @@ void InputControlConstrain(int control_type,double cm_min, double cm_max); bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); - void InputDepthAverageAtBase(int enum_type,int average_enum_type); + void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=ElementsEnum); void InputDuplicate(int original_enum,int new_enum); void InputScale(int enum_type,double scale_factor); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4930) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4931) @@ -14,5 +14,4 @@ #include #include "../objects.h" -#include "../../EnumDefinitions/EnumDefinitions.h" #include "../../shared/shared.h" #include "../../Container/Container.h" @@ -629,5 +628,5 @@ } else if (control_type==RheologyBEnum){ - inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBEnum); + inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum); Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; } @@ -870,5 +869,5 @@ GradjDrag(gradient); } - else if (control_type==RheologyBEnum){ + else if (control_type==RheologyB2dEnum){ GradjB(gradient); } @@ -966,5 +965,5 @@ adjointx_input=inputs->GetInput(AdjointxEnum); adjointy_input=inputs->GetInput(AdjointyEnum); - rheologyb_input=inputs->GetInput(RheologyBEnum); + rheologyb_input=inputs->GetInput(RheologyB2dEnum); /* Start looping on the number of gaussian points: */ @@ -1231,5 +1230,5 @@ /*Move input to Material if required (needed if control method) TO BE IMPROVED*/ - if (YEnum==RheologyBEnum){ + if (YEnum==RheologyBEnum || YEnum==RheologyB2dEnum){ this->matice->inputs->AddInput((Input*)yinput->copy()); } @@ -1287,5 +1286,5 @@ /*}}}*/ /*FUNCTION Tria::InputDepthAverageAtBase {{{1*/ -void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){ +void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ /*New input*/ @@ -1294,5 +1293,10 @@ /*copy input of enum_type*/ - oldinput=this->inputs->GetInput(enum_type); + if (object_enum==ElementsEnum) + oldinput=(Input*)this->inputs->GetInput(enum_type); + else if (object_enum==MaterialsEnum) + oldinput=(Input*)this->matice->inputs->GetInput(enum_type); + else + ISSMERROR("object %s not supported yet",EnumAsString(object_enum)); if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumAsString(enum_type)); newinput=(Input*)oldinput->copy(); @@ -1302,5 +1306,10 @@ /*Add new input to current element*/ - this->inputs->AddInput(newinput); + if (object_enum==ElementsEnum) + this->inputs->AddInput((Input*)newinput); + else if (object_enum==MaterialsEnum) + this->matice->inputs->AddInput((Input*)newinput); + else + ISSMERROR("object %s not supported yet",EnumAsString(object_enum)); } @@ -2126,5 +2135,5 @@ if (iomodel->rheology_B) { for(i=0;i<3;i++)nodeinputs[i]=iomodel->rheology_B[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(RheologyBEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(RheologyB2dEnum,nodeinputs)); } if (iomodel->control_parameter) { @@ -6048,4 +6057,5 @@ name==VyEnum || name==RheologyBEnum || + name==RheologyB2dEnum || name==RheologyNEnum || name==FitEnum || Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 4930) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 4931) @@ -20,4 +20,5 @@ #include "../../include/include.h" #include "../../shared/Exceptions/exceptions.h" +#include "../../EnumDefinitions/EnumDefinitions.h" /*}}}*/ @@ -84,5 +85,5 @@ void InputControlConstrain(int control_type,double cm_min, double cm_max); bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums); - void InputDepthAverageAtBase(int enum_type,int average_enum_type); + void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=ElementsEnum); void InputDuplicate(int original_enum,int new_enum); void InputScale(int enum_type,double scale_factor); Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 4930) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 4931) @@ -23,5 +23,5 @@ } /*}}}*/ -/*FUNCTION DoubleInput::DoubleInput(double* values){{{1*/ +/*FUNCTION DoubleInput::DoubleInput(double value){{{1*/ DoubleInput::DoubleInput(int in_enum_type,IssmDouble in_value){ @@ -253,2 +253,45 @@ } /*}}}*/ +/*FUNCTION DoubleInput::VerticallyIntegrate{{{1*/ +void DoubleInput::VerticallyIntegrate(Input* thickness_input){ + + /*Intermediaries*/ + double thickness_value; + + /*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())); + + /*vertically integrate depending on type:*/ + switch(thickness_input->Enum()){ + + case PentaVertexInputEnum: + thickness_input->GetParameterAverage(&thickness_value); + this->value=this->value*thickness_value; + return; + + default: + ISSMERROR("not implemented yet"); + } +} +/*}}}*/ +/*FUNCTION DoubleInput::PointwiseDivide{{{1*/ +Input* DoubleInput::PointwiseDivide(Input* inputB){ + + /*Ouput*/ + DoubleInput* outinput=NULL; + + /*Intermediaries*/ + double AdotBvalue; + double Bvalue; + + /*Check that inputB is of the same type*/ + inputB->GetParameterAverage(&Bvalue); + + /*Create new DoubleInput*/ + outinput=new DoubleInput(this->enum_type,this->value/Bvalue); + + /*Return output pointer*/ + return outinput; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 4930) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 4931) @@ -38,5 +38,5 @@ int EnumType(); Input* SpawnTriaInput(int* indices); - Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");}; + Input* PointwiseDivide(Input* inputB); ElementResult* SpawnResult(int step, double time); /*}}}*/ @@ -62,5 +62,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 VerticallyIntegrate(Input* thickness_input); void GetVectorFromInputs(Vec vector,int* doflist); void GetValuesPtr(double** pvalues,int* pnum_values); Index: /issm/trunk/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 4930) +++ /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 4931) @@ -47,5 +47,5 @@ if (iomodel->rheology_B) { for(i=0;irheology_B[int(iomodel->elements[num_vertices*index+i]-1)]; - this->inputs->AddInput(new TriaVertexInput(RheologyBEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(RheologyB2dEnum,nodeinputs)); } @@ -245,4 +245,14 @@ } /*}}}*/ +/*FUNCTION Matice::GetB2d {{{1*/ +double Matice::GetB2d(){ + + /*Output*/ + double B2d; + + inputs->GetParameterAverage(&B2d,RheologyB2dEnum); + return B2d; +} +/*}}}*/ /*FUNCTION Matice::GetN {{{1*/ double Matice::GetN(){ @@ -280,5 +290,5 @@ /*Get B and n*/ - B=GetB(); + B=GetB2d(); n=GetN(); @@ -488,5 +498,5 @@ /*Get B and n*/ - B=GetB(); + B=GetB2d(); n=GetN(); Index: /issm/trunk/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.h (revision 4930) +++ /issm/trunk/src/c/objects/Materials/Matice.h (revision 4931) @@ -58,4 +58,5 @@ void GetViscosityComplement(double* pviscosity_complement, double* pepsilon); double GetB(); + double GetB2d(); double GetN(); /*}}}*/