Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 15416) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 15417) @@ -533,6 +533,6 @@ ./classes/ElementResults/PentaP1ElementResult.h\ ./classes/ElementResults/PentaP1ElementResult.cpp\ - ./classes/Inputs/PentaP1Input.h\ - ./classes/Inputs/PentaP1Input.cpp\ + ./classes/Inputs/PentaInput.h\ + ./classes/Inputs/PentaInput.cpp\ ./classes/Elements/Penta.h\ ./classes/Elements/Penta.cpp\ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15417) @@ -182,5 +182,5 @@ if (!IsOnBed() || IsFloating()){ //empty friction: - this->inputs->AddInput(new PentaP1Input(BasalFrictionEnum,&basalfriction[0])); + this->inputs->AddInput(new PentaInput(BasalFrictionEnum,&basalfriction[0],P1Enum)); return; } @@ -209,5 +209,5 @@ /*Create PentaVertex input, which will hold the basal friction:*/ - this->inputs->AddInput(new PentaP1Input(BasalFrictionEnum,&basalfriction[0])); + this->inputs->AddInput(new PentaInput(BasalFrictionEnum,&basalfriction[0],P1Enum)); /*Clean up and return*/ @@ -354,10 +354,10 @@ /*Add Stress tensor components into inputs*/ - this->inputs->AddInput(new PentaP1Input(StressTensorxxEnum,&sigma_xx[0])); - this->inputs->AddInput(new PentaP1Input(StressTensorxyEnum,&sigma_xy[0])); - this->inputs->AddInput(new PentaP1Input(StressTensorxzEnum,&sigma_xz[0])); - this->inputs->AddInput(new PentaP1Input(StressTensoryyEnum,&sigma_yy[0])); - this->inputs->AddInput(new PentaP1Input(StressTensoryzEnum,&sigma_yz[0])); - this->inputs->AddInput(new PentaP1Input(StressTensorzzEnum,&sigma_zz[0])); + this->inputs->AddInput(new PentaInput(StressTensorxxEnum,&sigma_xx[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensorxyEnum,&sigma_xy[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensorxzEnum,&sigma_xz[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensoryyEnum,&sigma_yy[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum)); /*Clean up and return*/ @@ -763,9 +763,9 @@ for (int imonth=0;imonth<12;imonth++) { for(i=0;iAddTimeInput(newmonthinput1,time+imonth/12.*yts); for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); } @@ -1420,5 +1420,5 @@ /*create static input: */ - this->inputs->AddInput(new PentaP1Input(vector_enum,nodeinputs)); + this->inputs->AddInput(new PentaInput(vector_enum,nodeinputs,P1Enum)); } else if(M==numberofvertices+1){ @@ -1436,5 +1436,5 @@ if(t==0)transientinput=new TransientInput(vector_enum); - transientinput->AddTimeInput(new PentaP1Input(vector_enum,nodeinputs),time); + transientinput->AddTimeInput(new PentaInput(vector_enum,nodeinputs,P1Enum),time); } this->inputs->AddInput(transientinput); @@ -1490,5 +1490,5 @@ /*OK, we are on bed. Initialize global inputs as 0*/ - total_thickness_input =new PentaP1Input(ThicknessEnum,zeros_list); + total_thickness_input =new PentaInput(ThicknessEnum,zeros_list,P1Enum); /*Now follow all the upper element from the base to the surface to integrate the input*/ @@ -1508,8 +1508,8 @@ /*If first time, initialize total_integrated_input*/ if (step==0){ - if (original_input->ObjectEnum()==PentaP1InputEnum) - total_integrated_input=new PentaP1Input(average_enum_type,zeros_list); + if (original_input->ObjectEnum()==PentaInputEnum) + total_integrated_input=new PentaInput(average_enum_type,zeros_list,P1Enum); else if (original_input->ObjectEnum()==ControlInputEnum) - total_integrated_input=new PentaP1Input(average_enum_type,zeros_list); + total_integrated_input=new PentaInput(average_enum_type,zeros_list,P1Enum); else if (original_input->ObjectEnum()==DoubleInputEnum) total_integrated_input=new DoubleInput(average_enum_type,0.0); @@ -1524,5 +1524,5 @@ Helem_list[i+3]=Helem_list[i]; } - element_thickness_input=new PentaP1Input(ThicknessEnum,Helem_list); + element_thickness_input=new PentaInput(ThicknessEnum,Helem_list,P1Enum); /*Step3: Vertically integrate A COPY of the original*/ @@ -1752,5 +1752,5 @@ for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1760,5 +1760,5 @@ for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VxEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1768,5 +1768,5 @@ for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VyEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1776,5 +1776,5 @@ for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1830,5 +1830,5 @@ for(i=0;iData(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_cm_responses+i]; - datasetinput->inputs->AddObject(new PentaP1Input(InversionCostFunctionsCoefficientsEnum,nodeinputs)); + datasetinput->inputs->AddObject(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum)); } @@ -1983,7 +1983,7 @@ for(;;){ /*Add input to the element: */ - penta->inputs->AddInput(new PentaP1Input(ThicknessEnum,newthickness)); - penta->inputs->AddInput(new PentaP1Input(SurfaceEnum,newsurface)); - penta->inputs->AddInput(new PentaP1Input(BedEnum,newbed)); + penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum)); + penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum)); + penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum)); /*Stop if we have reached the surface*/ @@ -2016,5 +2016,5 @@ /*Add input to the element: */ - this->inputs->AddInput(new PentaP1Input(enum_type,values)); + this->inputs->AddInput(new PentaInput(enum_type,values,P1Enum)); /*Free ressources:*/ @@ -2049,5 +2049,5 @@ for(;;){ /*Add input to the element: */ - penta->inputs->AddInput(new PentaP1Input(enum_type,values)); + penta->inputs->AddInput(new PentaInput(enum_type,values,P1Enum)); /*Stop if we have reached the surface*/ @@ -2079,8 +2079,8 @@ /*update input*/ if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new PentaP1Input(name,values)); + material->inputs->AddInput(new PentaInput(name,values,P1Enum)); } else{ - this->inputs->AddInput(new PentaP1Input(name,values)); + this->inputs->AddInput(new PentaInput(name,values,P1Enum)); } return; @@ -2092,8 +2092,8 @@ /*update input*/ if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new PentaP1Input(name,values)); + material->inputs->AddInput(new PentaInput(name,values,P1Enum)); } else{ - this->inputs->AddInput(new PentaP1Input(name,values)); + this->inputs->AddInput(new PentaInput(name,values,P1Enum)); } return; @@ -2109,5 +2109,5 @@ } /*Add input to the element: */ - this->inputs->AddInput(new PentaP1Input(name,values)); + this->inputs->AddInput(new PentaInput(name,values,P1Enum)); /*Free ressources:*/ @@ -2121,5 +2121,5 @@ } /*Add input to the element: */ - this->inputs->AddInput(new PentaP1Input(name,values)); + this->inputs->AddInput(new PentaInput(name,values,P1Enum)); /*Free ressources:*/ @@ -2464,5 +2464,5 @@ /*Update inputs*/ - this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); + this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum); @@ -2800,5 +2800,5 @@ } //end of the loop over the vertices /*Update inputs*/ - this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0])); + this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&smb[0],P1Enum)); } /*}}}*/ @@ -2974,26 +2974,26 @@ if(!iomodel->Data(VxEnum)){ for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaP1Input(VxEnum,nodeinputs)); - if(dakota_analysis) this->inputs->AddInput(new PentaP1Input(QmuVxEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VxEnum,nodeinputs,P1Enum)); + if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVxEnum,nodeinputs,P1Enum)); } if(!iomodel->Data(VyEnum)){ for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaP1Input(VyEnum,nodeinputs)); - if(dakota_analysis) this->inputs->AddInput(new PentaP1Input(QmuVyEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VyEnum,nodeinputs,P1Enum)); + if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVyEnum,nodeinputs,P1Enum)); } if(!iomodel->Data(VzEnum)){ for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaP1Input(VzEnum,nodeinputs)); - if(dakota_analysis) this->inputs->AddInput(new PentaP1Input(QmuVzEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VzEnum,nodeinputs,P1Enum)); + if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVzEnum,nodeinputs,P1Enum)); } if(!iomodel->Data(PressureEnum)){ for(i=0;i<6;i++)nodeinputs[i]=0; if(dakota_analysis){ - this->inputs->AddInput(new PentaP1Input(PressureEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(QmuPressureEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(QmuPressureEnum,nodeinputs,P1Enum)); } if(isstokes){ - this->inputs->AddInput(new PentaP1Input(PressureEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(PressurePicardEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(PressurePicardEnum,nodeinputs,P1Enum)); } } @@ -3002,12 +3002,12 @@ if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){ for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum)); for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]); - this->inputs->AddInput(new PentaP1Input(VzPattynEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum)); } else{ for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(VzPattynEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum)); } } @@ -3016,12 +3016,12 @@ if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){ for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum)); for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]); - this->inputs->AddInput(new PentaP1Input(VzMacAyealEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum)); } else{ for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(VzMacAyealEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum)); } } @@ -3031,11 +3031,11 @@ /*Initialize mesh velocity*/ for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaP1Input(VxMeshEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(VyMeshEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(VzMeshEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VxMeshEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(VyMeshEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(VzMeshEnum,nodeinputs,P1Enum)); if(dakota_analysis){ - this->inputs->AddInput(new PentaP1Input(QmuVxMeshEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(QmuVyMeshEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(QmuVzMeshEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(QmuVxMeshEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(QmuVyMeshEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(QmuVzMeshEnum,nodeinputs,P1Enum)); } break; @@ -3044,7 +3044,7 @@ /*Initialize mesh velocity*/ for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaP1Input(VxMeshEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(VyMeshEnum,nodeinputs)); - this->inputs->AddInput(new PentaP1Input(VzMeshEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(VxMeshEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(VyMeshEnum,nodeinputs,P1Enum)); + this->inputs->AddInput(new PentaInput(VzMeshEnum,nodeinputs,P1Enum)); if (iomodel->Data(TemperatureEnum) && iomodel->Data(WaterfractionEnum) && iomodel->Data(PressureEnum)) { for(i=0;i<6;i++){ @@ -3056,5 +3056,5 @@ +latentheat*iomodel->Data(WaterfractionEnum)[penta_vertex_ids[i]-1]; } - this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(EnthalpyEnum,nodeinputs,P1Enum)); } else _error_("temperature and waterfraction required for the enthalpy solution"); @@ -3114,5 +3114,5 @@ /*Create PentaVertex input, which will hold the basal friction:*/ - this->inputs->AddInput(new PentaP1Input(ViscousHeatingEnum,&viscousheating[0])); + this->inputs->AddInput(new PentaInput(ViscousHeatingEnum,&viscousheating[0],P1Enum)); /*Clean up and return*/ @@ -4399,5 +4399,5 @@ this->inputs->GetInputValue(&converged,ConvergedEnum); if(converged){ - this->inputs->AddInput(new PentaP1Input(TemperatureEnum,values)); + this->inputs->AddInput(new PentaInput(TemperatureEnum,values,P1Enum)); /*Update Rheology only if converged (we must make sure that the temperature is below melting point @@ -4411,5 +4411,5 @@ B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0); for(i=0;imaterial->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); break; case ArrheniusEnum: @@ -4419,5 +4419,5 @@ material->GetN()); for(i=0;imaterial->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); break; default: @@ -4427,5 +4427,5 @@ } else{ - this->inputs->AddInput(new PentaP1Input(TemperaturePicardEnum,values)); + this->inputs->AddInput(new PentaInput(TemperaturePicardEnum,values,P1Enum)); } @@ -4475,7 +4475,7 @@ } - this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,values)); - this->inputs->AddInput(new PentaP1Input(WaterfractionEnum,waterfraction)); - this->inputs->AddInput(new PentaP1Input(TemperatureEnum,temperatures)); + this->inputs->AddInput(new PentaInput(EnthalpyEnum,values,P1Enum)); + this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum)); + this->inputs->AddInput(new PentaInput(TemperatureEnum,temperatures,P1Enum)); /*Update Rheology only if converged (we must make sure that the temperature is below melting point @@ -4489,5 +4489,5 @@ B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0); for(i=0;imaterial->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); break; case ArrheniusEnum: @@ -4497,5 +4497,5 @@ material->GetN()); for(i=0;imaterial->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); break; default: @@ -4505,5 +4505,5 @@ } else{ - this->inputs->AddInput(new PentaP1Input(EnthalpyPicardEnum,values)); + this->inputs->AddInput(new PentaInput(EnthalpyPicardEnum,values,P1Enum)); } @@ -4581,5 +4581,5 @@ GradientIndexing(&vertexpidlist[0],control_index); for(int i=0;iSetGradient(grad_input); @@ -5255,8 +5255,8 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(AdjointxEnum,lambdax)); - this->inputs->AddInput(new PentaP1Input(AdjointyEnum,lambday)); - this->inputs->AddInput(new PentaP1Input(AdjointzEnum,lambdaz)); - this->inputs->AddInput(new PentaP1Input(AdjointpEnum,lambdap)); + this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum)); + this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum)); + this->inputs->AddInput(new PentaInput(AdjointzEnum,lambdaz,P1Enum)); + this->inputs->AddInput(new PentaInput(AdjointpEnum,lambdap,P1Enum)); /*Free ressources:*/ @@ -5292,6 +5292,6 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(AdjointxEnum,lambdax)); - this->inputs->AddInput(new PentaP1Input(AdjointyEnum,lambday)); + this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum)); + this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum)); /*Free ressources:*/ @@ -5583,5 +5583,5 @@ values[i]=vector[vertexpidlist[i]]; } - new_input = new PentaP1Input(control_enum,values); + new_input = new PentaInput(control_enum,values,P1Enum); if(control_enum==MaterialsRheologyBbarEnum){ @@ -5614,5 +5614,5 @@ case VertexEnum: - /*New PentaP1Input*/ + /*New PentaInput*/ IssmDouble values[6]; @@ -5687,12 +5687,12 @@ /*Add new inputs: */ - this->inputs->AddInput(new PentaP1Input(ThicknessEnum,thickness)); - this->inputs->AddInput(new PentaP1Input(BedEnum,bed)); - this->inputs->AddInput(new PentaP1Input(SurfaceEnum,surface)); + this->inputs->AddInput(new PentaInput(ThicknessEnum,thickness,P1Enum)); + this->inputs->AddInput(new PentaInput(BedEnum,bed,P1Enum)); + this->inputs->AddInput(new PentaInput(SurfaceEnum,surface,P1Enum)); /*}}}*/ break; default: - this->inputs->AddInput(new PentaP1Input(name,values)); + this->inputs->AddInput(new PentaInput(name,values,P1Enum)); } break; @@ -5747,5 +5747,5 @@ if(t==0) transientinput=new TransientInput(name); - transientinput->AddTimeInput(new PentaP1Input(name,values),time); + transientinput->AddTimeInput(new PentaInput(name,values,P1Enum),time); transientinput->Configure(parameters); } @@ -8439,8 +8439,8 @@ /*Add vx and vy as inputs to the tria element: */ - penta->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - penta->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - penta->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - penta->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Stop if we have reached the surface*/ @@ -8529,8 +8529,8 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - this->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - this->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Free ressources:*/ @@ -8604,5 +8604,5 @@ Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum); if (vzmacayeal_input){ - if (vzmacayeal_input->ObjectEnum()!=PentaP1InputEnum){ + if (vzmacayeal_input->ObjectEnum()!=PentaInputEnum){ _error_("Cannot compute Vel as VzMacAyeal is of type " << EnumToStringx(vzmacayeal_input->ObjectEnum())); } @@ -8627,10 +8627,10 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - this->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - this->inputs->AddInput(new PentaP1Input(VzEnum,vz)); - this->inputs->AddInput(new PentaP1Input(VzStokesEnum,vzstokes)); - this->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum)); + this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum)); + this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Free ressources:*/ @@ -8704,8 +8704,8 @@ /*Add vx and vy as inputs to the tria element: */ - penta->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - penta->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - penta->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - penta->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Stop if we have reached the surface*/ @@ -8785,8 +8785,8 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - this->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - this->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Free ressources:*/ @@ -8853,5 +8853,5 @@ Input* vzpattyn_input=inputs->GetInput(VzPattynEnum); if (vzpattyn_input){ - if (vzpattyn_input->ObjectEnum()!=PentaP1InputEnum){ + if (vzpattyn_input->ObjectEnum()!=PentaInputEnum){ _error_("Cannot compute Vel as VzPattyn is of type " << EnumToStringx(vzpattyn_input->ObjectEnum())); } @@ -8876,10 +8876,10 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - this->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - this->inputs->AddInput(new PentaP1Input(VzEnum,vz)); - this->inputs->AddInput(new PentaP1Input(VzStokesEnum,vzstokes)); - this->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum)); + this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum)); + this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Free ressources:*/ @@ -8942,8 +8942,8 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - this->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - this->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Free ressources:*/ @@ -8999,5 +8999,5 @@ Input* vzstokes_input=inputs->GetInput(VzStokesEnum); if (vzstokes_input){ - if (vzstokes_input->ObjectEnum()!=PentaP1InputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum())); + if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum())); GetInputListOnVertices(&vzstokes[0],VzStokesEnum); } @@ -9011,5 +9011,5 @@ Input* vzstokes_input=inputs->GetInput(VzStokesEnum); if (vzstokes_input){ - if (vzstokes_input->ObjectEnum()!=PentaP1InputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum())); + if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum())); GetInputListOnVertices(&vzstokes[0],VzStokesEnum); } @@ -9039,14 +9039,14 @@ if(approximation!=PattynStokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum){ this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum); - this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); } else if(approximation==PattynStokesApproximationEnum){ - this->inputs->AddInput(new PentaP1Input(VzPattynEnum,vzpattyn)); + this->inputs->AddInput(new PentaInput(VzPattynEnum,vzpattyn,P1Enum)); } else if(approximation==MacAyealStokesApproximationEnum){ - this->inputs->AddInput(new PentaP1Input(VzMacAyealEnum,vzmacayeal)); - } - this->inputs->AddInput(new PentaP1Input(VzEnum,vz)); - this->inputs->AddInput(new PentaP1Input(VelEnum,vel)); + this->inputs->AddInput(new PentaInput(VzMacAyealEnum,vzmacayeal,P1Enum)); + } + this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum)); + this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); /*Free ressources:*/ @@ -9105,9 +9105,9 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaP1Input(VxEnum,vx)); - this->inputs->AddInput(new PentaP1Input(VyEnum,vy)); - this->inputs->AddInput(new PentaP1Input(VzEnum,vz)); - this->inputs->AddInput(new PentaP1Input(VelEnum,vel)); - this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum)); + this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum)); + this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum)); + this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum)); + this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum)); /*Free ressources:*/ @@ -9350,6 +9350,6 @@ for(;;){ /*Add input to the element: */ - penta->inputs->AddInput(new PentaP1Input(SedimentHeadEnum,values)); - penta->inputs->AddInput(new PentaP1Input(SedimentHeadResidualEnum,residual)); + penta->inputs->AddInput(new PentaInput(SedimentHeadEnum,values,P1Enum)); + penta->inputs->AddInput(new PentaInput(SedimentHeadResidualEnum,residual,P1Enum)); /*Stop if we have reached the surface*/ @@ -9451,10 +9451,10 @@ if(!this->IsFloating() && floatingelement==true){ for(i=0;iinputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0])); + this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum,&melting[0],P1Enum)); } /*Update inputs*/ - this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0])); - this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0])); + this->inputs->AddInput(new PentaInput(SurfaceEnum,&s[0],P1Enum)); + this->inputs->AddInput(new PentaInput(BedEnum,&b[0],P1Enum)); this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement)); @@ -9462,5 +9462,5 @@ if(migration_style==SubelementMigrationEnum){ for(i=0;iinputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0])); + this->inputs->AddInput(new PentaInput(GLlevelsetEnum,&phi[0],P1Enum)); this->InputExtrude(GLlevelsetEnum,ElementEnum); } Index: /issm/trunk-jpl/src/c/classes/Inputs/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/CMakeLists.txt (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Inputs/CMakeLists.txt (revision 15417) @@ -13,4 +13,4 @@ # }}} # THREED_SOURCES {{{ -set(THREED_SOURCES $ENV{ISSM_DIR}/src/c/classes/objects/Inputs/PentaP1Input.cpp PARENT_SCOPE) +set(THREED_SOURCES $ENV{ISSM_DIR}/src/c/classes/objects/Inputs/PentaInput.cpp PARENT_SCOPE) # }}} Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 15417) @@ -36,9 +36,9 @@ maxvalues =new TriaInput(enum_type,pmax,P1Enum); break; - case PentaP1InputEnum: - values =new PentaP1Input(enum_type,pvalues); - savedvalues=new PentaP1Input(enum_type,pvalues); - minvalues =new PentaP1Input(enum_type,pmin); - maxvalues =new PentaP1Input(enum_type,pmax); + case PentaInputEnum: + values =new PentaInput(enum_type,pvalues,P1Enum); + savedvalues=new PentaInput(enum_type,pvalues,P1Enum); + minvalues =new PentaInput(enum_type,pmin,P1Enum); + maxvalues =new PentaInput(enum_type,pmax,P1Enum); break; default: @@ -115,4 +115,8 @@ /*Object functions*/ +/*FUNCTION ControlInput::AXPY(){{{*/ +void ControlInput::AXPY(Input* xinput,IssmDouble scalar){ + values->AXPY(xinput,scalar); +}/*}}}*/ /*FUNCTION ControlInput::Constrain(){{{*/ void ControlInput::Constrain(void){ Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 15417) @@ -74,5 +74,5 @@ void Scale(IssmDouble scale_factor){_error_("not implemented yet");}; void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");}; - void AXPY(Input* xinput,IssmDouble scalar){_error_("not implemented yet");}; + void AXPY(Input* xinput,IssmDouble scalar); void Constrain(void); void Constrain(IssmDouble min,IssmDouble max); Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp (revision 15417) @@ -263,5 +263,5 @@ switch(thickness_input->ObjectEnum()){ - case PentaP1InputEnum: + case PentaInputEnum: thickness_input->GetInputAverage(&thickness_value); this->value=this->value*thickness_value; Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp (revision 15417) +++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp (revision 15417) @@ -0,0 +1,616 @@ +/*!\file PentaInput.c + * \brief: implementation of the PentaInput object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../classes.h" +#include "../../shared/shared.h" + +/*PentaInput constructors and destructor*/ +/*FUNCTION PentaInput::PentaInput(){{{*/ +PentaInput::PentaInput(){ + values = NULL; +} +/*}}}*/ +/*FUNCTION PentaInput::PentaInput(int in_enum_type,IssmDouble* values,int element_type){{{*/ +PentaInput::PentaInput(int in_enum_type,IssmDouble* in_values,int element_type_in) + :PentaRef(1) +{ + + /*Set PentaRef*/ + this->SetElementType(element_type_in,0); + this->element_type=element_type_in; + + /*Set Enum*/ + enum_type=in_enum_type; + + /*Set values*/ + this->values=xNew(this->NumberofNodes()); + for(int i=0;iNumberofNodes();i++) values[i]=in_values[i]; +} +/*}}}*/ +/*FUNCTION PentaInput::~PentaInput(){{{*/ +PentaInput::~PentaInput(){ + xDelete(this->values); +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION PentaInput::Echo {{{*/ +void PentaInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION PentaInput::DeepEcho{{{*/ +void PentaInput::DeepEcho(void){ + + _printf_("PentaInput:\n"); + _printf_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n"); + _printf_(" values: ["); + for(int i=0;iNumberofNodes();i++) _printf_(" "<values[i]); + _printf_("]\n"); +} +/*}}}*/ +/*FUNCTION PentaInput::Id{{{*/ +int PentaInput::Id(void){ return -1; } +/*}}}*/ +/*FUNCTION PentaInput::ObjectEnum{{{*/ +int PentaInput::ObjectEnum(void){ + + return PentaInputEnum; + +} +/*}}}*/ +/*FUNCTION PentaInput::copy{{{*/ +Object* PentaInput::copy() { + + return new PentaInput(this->enum_type,this->values,this->element_type); + +} +/*}}}*/ + +/*PentaInput management*/ +/*FUNCTION PentaInput::InstanceEnum{{{*/ +int PentaInput::InstanceEnum(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION PentaInput::SpawnTriaInput{{{*/ +Input* PentaInput::SpawnTriaInput(int* indices){ + + /*output*/ + TriaInput* outinput=NULL; + IssmDouble newvalues[3]; //Assume P1 interpolation only for now + + /*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 PentaInput::SpawnResult{{{*/ +ElementResult* PentaInput::SpawnResult(int step, IssmDouble time){ + + return new PentaP1ElementResult(this->enum_type,this->values,step,time); + +} +/*}}}*/ + +/*Object functions*/ +/*FUNCTION PentaInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/ +void PentaInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){ + + /*Call PentaRef function*/ + PentaRef::GetInputValue(pvalue,&values[0],gauss); + +} +/*}}}*/ +/*FUNCTION PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){{{*/ +void PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){ + + /*Call PentaRef function*/ + PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss); +} +/*}}}*/ +/*FUNCTION PentaInput::GetVxStrainRate3d{{{*/ +void PentaInput::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]; + + _assert_(this->NumberofNodes()==6); //Check Tria too + + /*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 PentaInput::GetVyStrainRate3d{{{*/ +void PentaInput::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]; + + _assert_(this->NumberofNodes()==6); //Check Tria too + + /*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 PentaInput::GetVzStrainRate3d{{{*/ +void PentaInput::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); + + _assert_(this->NumberofNodes()==6); //Check Tria too + + /*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 PentaInput::GetVxStrainRate3dPattyn{{{*/ +void PentaInput::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); + + _assert_(this->NumberofNodes()==6); //Check Tria too + + + /*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 PentaInput::GetVyStrainRate3dPattyn{{{*/ +void PentaInput::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); + + _assert_(this->NumberofNodes()==6); //Check Tria too + + + /*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 PentaInput::ChangeEnum{{{*/ +void PentaInput::ChangeEnum(int newenumtype){ + this->enum_type=newenumtype; +} +/*}}}*/ +/*FUNCTION PentaInput::GetInputAverage{{{*/ +void PentaInput::GetInputAverage(IssmDouble* pvalue){ + + int numnodes = this->NumberofNodes(); + IssmDouble numnodesd = reCast(numnodes); + IssmDouble value = 0.; + + for(int i=0;iNumberofNodes(); + IssmDouble squaremin; + + /*Now, figure out minimum of valuescopy: */ + squaremin=pow(this->values[0],2); + for(int i=1;ivalues[i],2)values[i],2); + } + /*Assign output pointers:*/ + *psquaremin=squaremin; +} +/*}}}*/ +/*FUNCTION PentaInput::ConstrainMin{{{*/ +void PentaInput::ConstrainMin(IssmDouble minimum){ + + int numnodes = this->NumberofNodes(); + for(int i=0;iNumberofNodes(); + + for(int i=0;inorm) norm=fabs(values[i]); + return norm; +} +/*}}}*/ +/*FUNCTION PentaInput::Max{{{*/ +IssmDouble PentaInput::Max(void){ + + int numnodes=this->NumberofNodes(); + IssmDouble max=values[0]; + + for(int i=1;imax) max=values[i]; + } + return max; +} +/*}}}*/ +/*FUNCTION PentaInput::MaxAbs{{{*/ +IssmDouble PentaInput::MaxAbs(void){ + + int numnodes=this->NumberofNodes(); + IssmDouble max=fabs(values[0]); + + for(int i=1;imax) max=fabs(values[i]); + } + return max; +} +/*}}}*/ +/*FUNCTION PentaInput::Min{{{*/ +IssmDouble PentaInput::Min(void){ + + const int numnodes=this->NumberofNodes(); + IssmDouble min=values[0]; + + for(int i=1;iNumberofNodes(); + IssmDouble min=fabs(values[0]); + + for(int i=1;iNumberofNodes(); + for(int i=0;iNumberofNodes(); + PentaInput* xpentainput=NULL; + + /*If xinput is a ControlInput, take its values directly*/ + if(xinput->ObjectEnum()==ControlInputEnum){ + xinput=((ControlInput*)xinput)->values; + } + + /*xinput is of the same type, so cast it: */ + if(xinput->ObjectEnum()!=PentaInputEnum) + _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); + xpentainput=(PentaInput*)xinput; + if(xpentainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xpentainput->element_type)); + + /*Carry out the AXPY operation depending on type:*/ + for(int i=0;ivalues[i]=this->values[i]+scalar*xpentainput->values[i]; + +} +/*}}}*/ +/*FUNCTION PentaInput::Constrain{{{*/ +void PentaInput::Constrain(IssmDouble cm_min, IssmDouble cm_max){ + + int i; + const int numnodes=this->NumberofNodes(); + + 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 PentaInput::Extrude{{{*/ +void PentaInput::Extrude(void){ + + switch(this->element_type){ + case P1Enum: + for(int i=0;i<3;i++) this->values[3+i]=this->values[i]; + break; + default: + _error_("not supported yet for type "<element_type)); + } +} +/*}}}*/ +/*FUNCTION PentaInput::VerticallyIntegrate{{{*/ +void PentaInput::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(this->element_type){ + case P1Enum:{ + 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 supported yet for type "<element_type)); + } +} +/*}}}*/ +/*FUNCTION PentaInput::PointwiseDivide{{{*/ +Input* PentaInput::PointwiseDivide(Input* inputB){ + + /*Ouput*/ + PentaInput* outinput=NULL; + + /*Intermediaries*/ + PentaInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=PentaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(PentaInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); + + /*Allocate intermediary*/ + IssmDouble* AdotBvalues=xNew(numnodes); + + /*Create point wise sum*/ + for(int i=0;ivalues[i]!=0); + AdotBvalues[i]=this->values[i]/xinputB->values[i]; + } + + /*Create new Penta vertex input (copy of current input)*/ + outinput=new PentaInput(this->enum_type,AdotBvalues,this->element_type); + + /*Return output pointer*/ + xDelete(AdotBvalues); + return outinput; + +} +/*}}}*/ +/*FUNCTION PentaInput::PointwiseMin{{{*/ +Input* PentaInput::PointwiseMin(Input* inputB){ + + /*Ouput*/ + PentaInput* outinput=NULL; + + /*Intermediaries*/ + int i; + PentaInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + IssmDouble *minvalues = xNew(numnodes); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=PentaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(PentaInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); + + /*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 PentaInput(this->enum_type,&minvalues[0],this->element_type); + + /*Return output pointer*/ + xDelete(minvalues); + return outinput; +} +/*}}}*/ +/*FUNCTION PentaInput::PointwiseMax{{{*/ +Input* PentaInput::PointwiseMax(Input* inputB){ + + /*Ouput*/ + PentaInput* outinput=NULL; + + /*Intermediaries*/ + int i; + PentaInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + IssmDouble *maxvalues = xNew(numnodes); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=PentaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(PentaInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); + + /*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 PentaInput(this->enum_type,&maxvalues[0],this->element_type); + + /*Return output pointer*/ + xDelete(maxvalues); + return outinput; +} +/*}}}*/ +/*FUNCTION PentaInput::GetVectorFromInputs{{{*/ +void PentaInput::GetVectorFromInputs(Vector* vector,int* doflist){ + + const int numnodes=this->NumberofNodes(); + vector->SetValues(numnodes,doflist,this->values,INS_VAL); + +} /*}}}*/ +/*FUNCTION PentaInput::Configure{{{*/ +void PentaInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h (revision 15417) +++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h (revision 15417) @@ -0,0 +1,84 @@ +/*! \file PentaInput.h + * \brief: header file for PentaInput object + */ + +#ifndef _PENTAINPUT_H_ +#define _PENTAINPUT_H_ + +/*Headers:*/ +/*{{{*/ +#include "./Input.h" +#include "../Elements/PentaRef.h" +class GaussTria; +class GaussPenta; +/*}}}*/ + +class PentaInput: public Input, public PentaRef{ + + public: + int enum_type; + IssmDouble* values; + + /*PentaInput constructors, destructors*/ + PentaInput(); + PentaInput(int enum_type,IssmDouble* values,int element_type_in); + ~PentaInput(); + + /*Object virtual functions definitions */ + void Echo(); + void DeepEcho(); + int Id(); + int ObjectEnum(); + Object *copy(); + + /*PentaInput management*/ + int InstanceEnum(); + Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB); + Input* PointwiseMin(Input* inputB); + Input* PointwiseMax(Input* inputB); + ElementResult* SpawnResult(int step, IssmDouble time); + void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; + void Configure(Parameters* parameters); + /*numerics*/ + void GetInputValue(bool* pvalue){_error_("not implemented yet");}; + void GetInputValue(int* pvalue){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); + void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");}; + void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; + void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss); + void GetInputAverage(IssmDouble* pvalue); + void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");}; + void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){_error_("not implemented yet");}; + + void GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; + void GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; + void GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss); + void GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss); + void GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss); + void GetVxStrainRate3dPattyn(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss); + void GetVyStrainRate3dPattyn(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss); + void ChangeEnum(int newenumtype); + + void SquareMin(IssmDouble* psquaremin,Parameters* parameters); + void ConstrainMin(IssmDouble minimum); + void Scale(IssmDouble scale_factor); + void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");}; + void AXPY(Input* xinput,IssmDouble scalar); + void Constrain(IssmDouble cm_min, IssmDouble cm_max); + IssmDouble InfinityNorm(void); + IssmDouble Max(void); + IssmDouble MaxAbs(void); + IssmDouble Min(void); + IssmDouble MinAbs(void); + void Extrude(void); + void VerticallyIntegrate(Input* thickness_input); + void GetVectorFromInputs(Vector* vector,int* doflist); + +}; +#endif /* _PENTAINPUT_H */ Index: sm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.cpp (revision 15416) +++ (revision ) @@ -1,600 +1,0 @@ -/*!\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 "../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;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(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;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(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;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(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;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(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;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(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;ivalues[i]; - - /*Now, figure out minimum of valuescopy: */ - squaremin=pow(valuescopy[0],2); - for(i=1;inorm) 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;imax) 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;imax) 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;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) _error_("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: - _error_("not implemented yet"); - } - -} -/*}}}*/ -/*FUNCTION PentaP1Input::Constrain{{{*/ -void PentaP1Input::Constrain(IssmDouble cm_min, IssmDouble 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){ - - 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;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; - 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;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; - 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;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::Configure{{{*/ -void PentaP1Input::Configure(Parameters* parameters){ - /*do nothing: */ -} -/*}}}*/ Index: sm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.h (revision 15416) +++ (revision ) @@ -1,87 +1,0 @@ -/*! \file PentaP1Input.h - * \brief: header file for PentaP1Input object - */ - -#ifndef _PENTAP1INPUT_H_ -#define _PENTAP1INPUT_H_ - -/*Headers:*/ -/*{{{*/ -#include "./Input.h" -#include "../Elements/PentaRef.h" -class GaussTria; -class GaussPenta; -/*}}}*/ - -class PentaP1Input: public Input, public PentaRef{ - - public: - /*just hold 6 values for 6 vertices: */ - int enum_type; - IssmDouble values[6]; - - /*PentaP1Input constructors, destructors: {{{*/ - PentaP1Input(); - PentaP1Input(int enum_type,IssmDouble* values); - ~PentaP1Input(); - /*}}}*/ - /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); - int ObjectEnum(); - Object *copy(); - /*}}}*/ - /*PentaP1Input management: {{{*/ - int InstanceEnum(); - Input* SpawnTriaInput(int* indices); - Input* PointwiseDivide(Input* inputB); - Input* PointwiseMin(Input* inputB); - Input* PointwiseMax(Input* inputB); - ElementResult* SpawnResult(int step, IssmDouble time); - void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; - void Configure(Parameters* parameters); - /*}}}*/ - /*numerics: {{{*/ - void GetInputValue(bool* pvalue){_error_("not implemented yet");}; - void GetInputValue(int* pvalue){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss); - void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");}; - void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; - void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss); - void GetInputAverage(IssmDouble* pvalue); - void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");}; - void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){_error_("not implemented yet");}; - - void GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; - void GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");}; - void GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss); - void GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss); - void GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss); - void GetVxStrainRate3dPattyn(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss); - void GetVyStrainRate3dPattyn(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss); - void ChangeEnum(int newenumtype); - - void SquareMin(IssmDouble* psquaremin,Parameters* parameters); - void ConstrainMin(IssmDouble minimum); - void Scale(IssmDouble scale_factor); - void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");}; - void AXPY(Input* xinput,IssmDouble scalar); - void Constrain(IssmDouble cm_min, IssmDouble cm_max); - IssmDouble InfinityNorm(void); - IssmDouble Max(void); - IssmDouble MaxAbs(void); - IssmDouble Min(void); - IssmDouble MinAbs(void); - void Extrude(void); - void VerticallyIntegrate(Input* thickness_input); - void GetVectorFromInputs(Vector* vector,int* doflist); - /*}}}*/ - -}; -#endif /* _PENTAP1INPUT_H */ Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 15417) @@ -38,5 +38,4 @@ TriaInput::~TriaInput(){ xDelete(this->values); - return; } /*}}}*/ @@ -252,5 +251,5 @@ void TriaInput::ConstrainMin(IssmDouble minimum){ - const int numnodes = this->NumberofNodes(); + int numnodes = this->NumberofNodes(); for(int i=0;iNumberofNodes(); + IssmDouble norm=0.; + int numnodes=this->NumberofNodes(); for(int i=0;inorm) norm=fabs(values[i]); @@ -270,5 +269,5 @@ IssmDouble TriaInput::Max(void){ - const int numnodes=this->NumberofNodes(); + int numnodes=this->NumberofNodes(); IssmDouble max=values[0]; @@ -282,5 +281,5 @@ IssmDouble TriaInput::MaxAbs(void){ - const int numnodes=this->NumberofNodes(); + int numnodes=this->NumberofNodes(); IssmDouble max=fabs(values[0]); @@ -319,5 +318,4 @@ const int numnodes=this->NumberofNodes(); - for(int i=0;iNumberofNodes(); - TriaInput* xtriavertexinput=NULL; + TriaInput* xtriainput=NULL; /*xinput is of the same type, so cast it: */ if(xinput->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); - xtriavertexinput=(TriaInput*)xinput; - if(xtriavertexinput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); + xtriainput=(TriaInput*)xinput; + if(xtriainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); /*Carry out the AXPY operation depending on type:*/ - for(i=0;ivalues[i]=this->values[i]+scalar*xtriavertexinput->values[i]; + for(int i=0;ivalues[i]=this->values[i]+scalar*xtriainput->values[i]; } @@ -387,7 +384,7 @@ /*Check that inputB is of the same type*/ - if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); xinputB=(TriaInput*)inputB; - if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); /*Create point wise min*/ @@ -421,5 +418,5 @@ if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); xinputB=(TriaInput*)inputB; - if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type)); /*Create point wise max*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 15417) @@ -80,5 +80,4 @@ void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");}; void GetVectorFromInputs(Vector* vector,int* doflist); - /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15417) @@ -817,5 +817,5 @@ if (iomodel->Data(MaterialsRheologyBEnum)) { for(i=0;iData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,nodeinputs,P1Enum)); } @@ -823,5 +823,5 @@ if (iomodel->Data(MaterialsRheologyNEnum)) { for(i=0;iData(MaterialsRheologyNEnum)[index]; - this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum)); } @@ -829,5 +829,5 @@ if (iomodel->Data(MaterialsRheologyZEnum)) { for(i=0;iData(MaterialsRheologyZEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(MaterialsRheologyZEnum,nodeinputs,P1Enum)); } @@ -843,5 +843,5 @@ for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -852,5 +852,5 @@ for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 15417) @@ -14,5 +14,5 @@ #include "../Inputs/Inputs.h" #include "../Inputs/TriaInput.h" -#include "../Inputs/PentaP1Input.h" +#include "../Inputs/PentaInput.h" #include "../Inputs/ControlInput.h" #include "../Elements/Element.h" @@ -571,5 +571,5 @@ IssmDouble valuesp[6]; for (int i=0;i<6;i++) valuesp[i]=vector[((Penta*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector - this->inputs->AddInput(new PentaP1Input(name,valuesp)); + this->inputs->AddInput(new PentaInput(name,valuesp,P1Enum)); return; } @@ -637,5 +637,5 @@ IssmDouble valuesp[6]; for (int i=0;i<6;i++) valuesp[i]=vector[((Penta*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector - this->inputs->AddInput(new PentaP1Input(name,valuesp)); + this->inputs->AddInput(new PentaInput(name,valuesp,P1Enum)); return; } @@ -751,5 +751,5 @@ if (iomodel->Data(MaterialsRheologyBEnum)) { for(i=0;iData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,nodeinputs,P1Enum)); } @@ -757,5 +757,5 @@ if (iomodel->Data(MaterialsRheologyNEnum)) { for(i=0;iData(MaterialsRheologyNEnum)[index]; - this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs)); + this->inputs->AddInput(new PentaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum)); } @@ -771,5 +771,5 @@ for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 15416) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 15417) @@ -56,5 +56,5 @@ #include "./Inputs/DoubleInput.h" #include "./Inputs/IntInput.h" -#include "./Inputs/PentaP1Input.h" +#include "./Inputs/PentaInput.h" #include "./Inputs/TriaInput.h" #include "./Inputs/ControlInput.h" Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15416) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15417) @@ -370,5 +370,5 @@ PenpairEnum, PentaEnum, - PentaP1InputEnum, + PentaInputEnum, ProfilerEnum, MatrixParamEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15417) @@ -370,5 +370,5 @@ case PenpairEnum : return "Penpair"; case PentaEnum : return "Penta"; - case PentaP1InputEnum : return "PentaP1Input"; + case PentaInputEnum : return "PentaInput"; case ProfilerEnum : return "Profiler"; case MatrixParamEnum : return "MatrixParam"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15416) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15417) @@ -376,5 +376,5 @@ else if (strcmp(name,"Penpair")==0) return PenpairEnum; else if (strcmp(name,"Penta")==0) return PentaEnum; - else if (strcmp(name,"PentaP1Input")==0) return PentaP1InputEnum; + else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; else if (strcmp(name,"Profiler")==0) return ProfilerEnum; else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;