Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23543) +++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23544) @@ -74,7 +74,8 @@ int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum}; int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum}; int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum}; - femmodel->InitTransientOutputx(&stackedinput[0],4); + //femmodel->InitTransientOutputx(&stackedinput[0],4); + femmodel->InitMeanOutputx(&stackedinput[0],4); while(hydrotimeStackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4); + //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4); + femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],4); } - femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4); + //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4); + femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],4); } else{ int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum}; int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum}; int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum}; - femmodel->InitTransientOutputx(&stackedinput[0],2); + //femmodel->InitTransientOutputx(&stackedinput[0],2); + femmodel->InitMeanOutputx(&stackedinput[0],2); while(hydrotimeStackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2); + //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2); + femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],2); } - femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2); + //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2); + femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],2); } } else{ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23543) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23544) @@ -4820,15 +4820,15 @@ if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n"); } /*}}}*/ -void FemModel::InitTransientOutputx(int* input_enum,int numoutputs){ /*{{{*/ +void FemModel::InitTransientOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/ for(int i=0;iSize();j++){ - TransientInput* transient_input = new TransientInput(input_enum[i]); + TransientInput* transient_input = new TransientInput(stackedinput_enum[i]); /*Intermediaries*/ Element* element=xDynamicCast(elements->GetObjectByOffset(j)); //transient_input->Configure(element->parameters); @@ -4849,13 +4849,12 @@ /*Intermediaries*/ Element* element=xDynamicCast(elements->GetObjectByOffset(j)); TransientInput* stacking_input=NULL; - Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input); - Input* input_to_stack=element->GetInput(input_enum[i]); _assert_(input_to_stack); + Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input); //this is the enum stack stacking_input=dynamic_cast(input); int numvertices = element->GetNumberOfVertices(); IssmDouble* N=xNew(numvertices); - element->GetInputListOnVertices(&N[0],input_enum[i]); + element->GetInputListOnVertices(&N[0],input_enum[i]); //this is the enum to stack switch(element->ObjectEnum()){ case TriaEnum: stacking_input->AddTimeInput(new TriaInput(stackedinput_enum[i],&N[0],P1Enum),hydrotime); @@ -4874,10 +4873,10 @@ } } /*}}}*/ -void FemModel::AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/ +void FemModel::AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/ for(int i=0;iGetNumberOfVertices(); IssmDouble* time_averaged=NULL; - Input* input=element->inputs->GetInput(input_enum[i]); _assert_(input); + Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input); TransientInput* stacking_input=NULL; stacking_input=dynamic_cast(input); stacking_input->GetInputAverageOnTimes(&time_averaged,init_time); @@ -4899,6 +4898,95 @@ } } /*}}}*/ +void FemModel::InitMeanOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/ + + for(int i=0;iSize();j++){ + /*Intermediaries*/ + Element* element =xDynamicCast(elements->GetObjectByOffset(j)); + int numvertices =element->GetNumberOfVertices(); + IssmDouble zeros[numvertices] ={0.0}; + switch(element->ObjectEnum()){ + case TriaEnum: + element->inputs->AddInput(new TriaInput(stackedinput_enum[i],&zeros[0],P1Enum)); + break; + case PentaEnum: + element->inputs->AddInput(new PentaInput(stackedinput_enum[i],&zeros[0],P1Enum)); + break; + case TetraEnum: + element->inputs->AddInput(new TetraInput(stackedinput_enum[i],&zeros[0],P1Enum)); + break; + default: _error_("Not implemented yet"); + } + } + } + } +} +/*}}}*/ +void FemModel::SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs){ /*{{{*/ + + //First get sub-timestep + IssmDouble hydrodt; + this->parameters->FindParam(&hydrodt,HydrologydtEnum); + + for(int i=0;iSize();j++){ + /*Intermediaries*/ + Element* element=xDynamicCast(elements->GetObjectByOffset(j)); + int numvertices = element->GetNumberOfVertices(); + IssmDouble* values_to_add=xNew(numvertices); + IssmDouble* existing_values=xNew(numvertices); + element->GetInputListOnVertices(&values_to_add[0],input_enum[i]); //those are the values to add + element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add + for(int k=0;kAddInput(stackedinput_enum[i],&existing_values[0],P1Enum); + xDelete(existing_values); + xDelete(values_to_add); + } + } + } +} +/*}}}*/ +void FemModel::AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs){ /*{{{*/ + + //First get timestep + IssmDouble maindt; + this->parameters->FindParam(&maindt,TimesteppingTimeStepEnum); + for(int i=0;iSize();j++){ + /*Intermediaries*/ + Element* element=xDynamicCast(elements->GetObjectByOffset(j)); + int numvertices = element->GetNumberOfVertices(); + IssmDouble* time_averaged=xNew(numvertices); + IssmDouble* existing_values=xNew(numvertices); + element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add + + for(int k=0;kAddInput(averagedinput_enum[i],&time_averaged[0],P1Enum); + xDelete(time_averaged); + xDelete(existing_values); + } + } + } +} +/*}}}*/ #ifdef _HAVE_JAVASCRIPT_ FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/ /*configuration: */ Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 23543) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 23544) @@ -140,14 +140,14 @@ void Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y); #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy); - void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy); + void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); #endif #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); + void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); IssmDouble SealevelriseOceanAverage(Vector* Sg); #endif void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); @@ -158,7 +158,10 @@ void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count); void InitTransientOutputx(int* input_enum, int numoutputs); void StackTransientOutputx(int* input_enum,int* stackedinput_enum, IssmDouble hydrotime, int numoutputs); - void AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble hydrotime,int numoutputs); + void AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs); + void InitMeanOutputx(int* stackedinput_enum,int numoutputs); + void SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs); + void AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs); void UpdateConstraintsx(void); int UpdateVertexPositionsx(void);