Changeset 25024 for issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

Timestamp:

06/12/20 08:22:16 (5 years ago)

Author:

bdef

Message:

CHG:moving averaging on inputs and modification on 905 test

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶

@@ -1009 +1009 @@
 /*}}}*/
 void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
-
         _assert_(end_time>start_time);
 
-        /*Intermediaries*/
-        IssmDouble averaged_values[NUMVERTICES];
-        IssmDouble current_values[NUMVERTICES];
-        IssmDouble dt;
-        int        found,start_offset,end_offset;
-
         /*Get transient input time steps*/
-        int         numtimesteps;
-        IssmDouble *timesteps    = NULL;
         TransientInput2* transient_input  = this->inputs2->GetTransientInput(transientinput_enum);
-        transient_input->GetAllTimes(&timesteps,&numtimesteps);
-
-        /*go through the timesteps, and grab offset for start and end*/
-        found=binary_search(&start_offset,start_time,timesteps,numtimesteps);
-        if(!found) _error_("Input not found (is TransientInput sorted ?)");
-        found=binary_search(&end_offset,end_time,timesteps,numtimesteps);
-        if(!found) _error_("Input not found (is TransientInput sorted ?)");
-
-        Gauss* gauss=this->NewGauss();
-
-        /*stack the input for each timestep in the slice*/
-        int offset = start_offset;
-        while(offset < end_offset ){
-
-
-                if(offset==-1){
-                        /*get values for the first time: */
-                        _assert_(start_time<timesteps[0]);
-                        PentaInput2* input = transient_input->GetPentaInput(0);
-
-                        int interpolation = input->GetInterpolation();
-                        _assert_(interpolation==P1Enum);
-                        /*Intermediaries*/
-                        int numindices;
-                        int indices[NUMVERTICES];
-                        numindices = NUMVERTICES;
-                        for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
-                        input->Serve(numindices,&indices[0]);
-
-                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                gauss->GaussVertex(iv);
-                                input->GetInputValue(&current_values[iv],gauss);
-                        }
-                        dt = timesteps[0] - start_time; _assert_(dt>0.);
-                }
-                else{
-                        PentaInput2* input = transient_input->GetPentaInput(offset+1);
-                        int interpolation = input->GetInterpolation();
-                        _assert_(interpolation==P1Enum);
-                        /*Intermediaries*/
-                        int numindices;
-                        int indices[NUMVERTICES];
-                        numindices = NUMVERTICES;
-                        for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
-                        input->Serve(numindices,&indices[0]);
-
-                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                gauss->GaussVertex(iv);
-                                input->GetInputValue(&current_values[iv],gauss);
-                        }
-                        if(offset == numtimesteps-1){
-                                dt = end_time - timesteps[offset]; _assert_(dt>0.);
-                        }
-                        else{
-                                dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
-                        }
-                }
-
-                switch(averaging_method){
-                        case 0: /*Arithmetic mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
-                                }
-                                break;
-                        case 1: /*Geometric mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] *= dt*current_values[iv];
-                                }
-                                break;
-                        case 2: /*Harmonic mean*/
-                                if(offset==start_offset){
-                                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                                _assert_(current_values[iv]>1.e-50);
-                                                averaged_values[iv]  = dt*1./current_values[iv];
-                                        }
-                                }
-                                else{
-                                        for(int iv=0;iv<NUMVERTICES;iv++){
-                                                _assert_(current_values[iv]>1.e-50);
-                                                averaged_values[iv]  += dt*1./current_values[iv];
-                                        }
-                                }
-                                break;
-                        default:
-                                _error_("averaging method is not recognised");
-                }
-
-                offset+=1;
-        }
-
-        /*Integration done, now normalize*/
-        switch(averaging_method){
-                case 0: //Arithmetic mean
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] =  averaged_values[iv]/(end_time - start_time);
-                        break;
-                case 1: /*Geometric mean*/
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = pow(averaged_values[iv], 1./(end_time - start_time));
-                        break;
-                case 2: /*Harmonic mean*/
-                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = 1./(averaged_values[iv]/(end_time - start_time));
-                        break;
-                default:
-                        _error_("averaging method is not recognised");
-        }
-
-        this->AddInput2(averagedinput_enum,&averaged_values[0],P1Enum);
-
-        /*Cleanup*/
-        delete gauss;
-        xDelete<IssmDouble>(timesteps);
+        PentaInput2* averaged_input = transient_input->GetPentaInput(start_time,end_time,averaging_method);
+        Input2* averaged_copy = averaged_input->copy();
+
+        averaged_input->ChangeEnum(averagedinput_enum);
+        this->inputs2->AddInput(averaged_copy);
 }
 /*}}}*/