Changeset 24366

Timestamp:

11/20/19 05:11:35 (5 years ago)

Author:

bdef

Message:

BUG: fixing temporal mean computation

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (14 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -280 +280 @@
         if(M!=iomodel->numberofvertices) _error_("not supported yet");
         if(N!=num_inputs) _error_("sizes are not consistent");
-        
+
 
         int        tria_vertex_ids[3];
@@ -2128 +2128 @@
         /*stack the input for each timestep in the slice*/
         int offset = start_offset;
-        while(offset <= end_offset ){
+        while(offset < end_offset ){
 
                 if(offset==-1){
@@ -2142 +2142 @@
                 }
                 else{
-                        TriaInput2* input = transient_input->GetTriaInput(offset);
+                        TriaInput2* input = transient_input->GetTriaInput(offset+1);
                         _assert_(input->GetInterpolation()==P1Enum);
                         this->InputServe(input);
@@ -2155 +2155 @@
                         dt = timesteps[0] - start_time;
                 }
-                else if(offset = numtimesteps-1){
+                else if(offset == numtimesteps-1){
                         dt = end_time - timesteps[offset];
                 }
@@ -2176 +2176 @@
                                 }
                                 else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
+                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] *= dt*current_values[iv];
                                 }
                                 break;
@@ -2189 +2189 @@
                                         for(int iv=0;iv<NUMVERTICES;iv++){
                                                 _assert_(current_values[iv]>1.e-50);
-                                                averaged_values[iv]  = dt*1./current_values[iv];
+                                                averaged_values[iv]  += dt*1./current_values[iv];
                                         }
                                 }
@@ -3843 +3843 @@
 
         /*First, serarch the input: */
-        Input2* data=this->GetInput2(natureofdataenum); 
+        Input2* data=this->GetInput2(natureofdataenum);
 
         /*figure out if we have the vertex id: */
@@ -4759 +4759 @@
         /*Get material parameters :*/
         rho_ice=FindParam(MaterialsRhoIceEnum);
-        Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input); 
+        Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input);
         Input2* gllevelset_input   = this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
         Input2* scalefactor_input  = NULL;
         if(scaled==true){
-                scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); 
+                scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
         }
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
@@ -4804 +4804 @@
         /*Get material parameters :*/
         rho_ice=FindParam(MaterialsRhoIceEnum);
-        Input2* groundedmelt_input = this->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input); 
+        Input2* groundedmelt_input = this->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input);
         Input2* gllevelset_input = this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
         Input2* scalefactor_input = NULL;
         if(scaled==true){
-                scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); 
+                scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
         }
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
@@ -5298 +5298 @@
 
         /*Compute ice thickness change: */
-        Input2* deltathickness_input=this->GetInput2(EsaDeltathicknessEnum); 
+        Input2* deltathickness_input=this->GetInput2(EsaDeltathicknessEnum);
         if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
         deltathickness_input->GetInputAverage(&I);
@@ -5436 +5436 @@
 
         /*Compute ice thickness change: */
-        Input2* deltathickness_input=this->GetInput2(EsaDeltathicknessEnum); 
+        Input2* deltathickness_input=this->GetInput2(EsaDeltathicknessEnum);
         if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
         deltathickness_input->GetInputAverage(&I);
@@ -5670 +5670 @@
 
                 /*Compute ice thickness change: */
-                Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum); 
+                Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum);
                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
                 deltathickness_input->GetInputAverage(&I);
@@ -5806 +5806 @@
 
         /*Compute ice thickness change: */
-        Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum); 
+        Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum);
         if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
 
@@ -6149 +6149 @@
 
         /*Compute ice thickness change: */
-        Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum); 
+        Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum);
         if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
         deltathickness_input->GetInputAverage(&I);