Changeset 22281

Timestamp:

11/30/17 08:06:43 (7 years ago)

Author:

bdef

Message:

CHG:implementing some 3D compliance and hydro_core clean-up

Location:

issm/trunk-jpl/src/c

Files:

• classes/Inputs/TransientInput.cpp (modified) (1 diff)
• cores/hydrology_core.cpp (modified) (6 diffs)

issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

@@ -315 +315 @@
 void TransientInput::GetInputAverageOnTimes(IssmDouble** pvalues, IssmDouble init_time){/*{{{*/
 
-        IssmDouble  time;
-        IssmDouble  preceding_time;
+        int                                     numnodes;
+        IssmDouble  subtimestep;
         IssmDouble* values= NULL;
-
-        values=xNew<IssmDouble>(3);
-        
+        Input* input=NULL;
+
+        /*Initialize numnode from transientinput out of time loop*/
         for(int i=0;i<this->numtimesteps;i++){
-                time                                             = this->timesteps[i];
-                TriaInput*      input            = (TriaInput*)this->inputs->GetObjectByOffset(i); _assert_(input);
-                int                                     numnodes = input->NumberofNodes(input->interpolation_type);
+                /*First compute the lengt of the sub-timestep*/
                 if(i==0){
-                        preceding_time = init_time;
-                        for(int j=0;j<numnodes;j++){
-                                values[j]=(input->values[j]*(time-preceding_time));
-                        }
+                        subtimestep = this->timesteps[i]-init_time;
                 }
                 else{
-                        preceding_time = this->timesteps[i-1];
+                        subtimestep = this->timesteps[i]-this->timesteps[i-1];
+                }       
+                /*Figure out type of input and get it*/
+                input = xDynamicCast<Input*>(this->inputs->GetObjectByOffset(i)); _assert_(input);
+                switch(input->ObjectEnum()){
+                case TriaInputEnum:{
+                        TriaInput* triainput = xDynamicCast<TriaInput*>(this->inputs->GetObjectByOffset(i)); _assert_(triainput);
+                        numnodes= triainput->NumberofNodes(triainput->interpolation_type);
+                        values=xNewZeroInit<IssmDouble>(numnodes);
+                        /*Sum the values weighted by their respective sub-timestep*/
                         for(int j=0;j<numnodes;j++){
-                                values[j]+=(input->values[j]*(time-preceding_time));
-                                if(i==this->numtimesteps-1){
-                                        values[j]=values[j]/(this->timesteps[this->numtimesteps-1]-init_time);
-                                }
+                                values[j]+=(triainput->values[j]*subtimestep);
                         }
+                        break;
                 }
+                case PentaInputEnum:{
+                        PentaInput*     pentainput = xDynamicCast<PentaInput*>(this->inputs->GetObjectByOffset(i)); _assert_(pentainput);
+                        numnodes= pentainput->NumberofNodes(pentainput->interpolation_type);
+                        values=xNewZeroInit<IssmDouble>(numnodes);
+                        /*Sum the values weighted by their respective sub-timestep*/
+                        for(int j=0;j<numnodes;j++){
+                                values[j]+=(pentainput->values[j]*subtimestep);
+                        }
+                        break;
+                }
+                case TetraInputEnum:{
+                        TetraInput*     tetrainput = xDynamicCast<TetraInput*>(this->inputs->GetObjectByOffset(i)); _assert_(tetrainput);
+                        numnodes= tetrainput->NumberofNodes(tetrainput->interpolation_type);
+                        values=xNewZeroInit<IssmDouble>(numnodes);
+                        /*Sum the values weighted by their respective sub-timestep*/
+                        for(int j=0;j<numnodes;j++){
+                                values[j]+=(tetrainput->values[j]*subtimestep);
+                        }
+                        break;
+                }
+                default: _error_("not implemented yet");
+                }
+        }
+        /*Compute the average based on the length of the full timestep*/
+        for(int j=0;j<numnodes;j++){
+                values[j]=values[j]/(this->timesteps[this->numtimesteps-1]-init_time);
         }
         *pvalues=values;

issm/trunk-jpl/src/c/cores/hydrology_core.cpp

@@ -24 +24 @@
         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);       
 
-        if(VerboseSolution()) _printf0_("   computing water heads\n");
-                        
-        /*first compute slopes: */
+        /*Using the Shreve based Model*/
         if (hydrology_model==HydrologyshreveEnum){
+                if(VerboseSolution()) _printf0_("   computing water heads\n");
+                /*first compute slopes: */
                 surfaceslope_core(femmodel);
                 bedslope_core(femmodel);
-        }
-
-        /*Using the Shreve based Model*/
-        if (hydrology_model==HydrologyshreveEnum){
+                /*and then go to water column*/
                 if(VerboseSolution()) _printf0_("   computing water column\n");
                 femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
@@ -45 +42 @@
                         int outputs[3] = {WatercolumnEnum,HydrologyWaterVxEnum,HydrologyWaterVyEnum};
                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
-                        
-                        /*unload results*/
-                        if(VerboseSolution()) _printf0_("   saving temporary results\n");
-                        OutputResultsx(femmodel);
                 }
         }
@@ -56 +49 @@
                 /*intermediary: */
                 int        step,hydroslices;
-                int        numoutputs;
+                int        numoutputs=0;
                 char       **requested_outputs = NULL;
                 IssmDouble time,init_time,hydrotime,yts;
@@ -111 +104 @@
                         if(VerboseSolution()) _printf0_("   saving results \n");
                         femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
-                        /*unload results removed 23/11 needs checking*/
-                        /* if(VerboseSolution()) _printf0_("   saving temporary results\n"); */
-                        /* OutputResultsx(femmodel); */
                 }
                 /*Free ressources:*/    
@@ -119 +109 @@
                         for(int i=0;i<numoutputs;i++){
                                 xDelete<char>(requested_outputs[i]);
-                        } 
+                        }
                         xDelete<char*>(requested_outputs);
                 }
@@ -138 +128 @@
                         int outputs[5] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum,HydrologyBasalFluxEnum,DegreeOfChannelizationEnum};
                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],5);
-                        
-                        /*unload results*/
-                        if(VerboseSolution()) _printf0_("   saving temporary results\n");
-                        OutputResultsx(femmodel);
                 }
         }