Changeset 21872


Ignore:
Timestamp:
07/26/17 11:15:15 (8 years ago)
Author:
Mathieu Morlighem
Message:

CHG: better way to depth average + some clean up

Location:
issm/trunk-jpl/src/c/classes
Files:
26 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r21811 r21872  
    216216                virtual int        GetNumberOfVertices(void)=0;
    217217                virtual void       GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int solutionenum)=0;
    218                 virtual Element*   GetUpperElement(void)=0;
    219218                virtual void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid)=0;
    220219                virtual void       GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;

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

    r21808 r21872  
    11271127}
    11281128/*}}}*/
    1129 Penta*     Penta::GetSurfacePenta(void){/*{{{*/
    1130 
    1131         /*Output*/
    1132         Penta* penta=NULL;
    1133 
    1134         /*Go through all pentas till the surface is reached*/
    1135         penta=this;
    1136         for(;;){
    1137                 /*Stop if we have reached the surface, else, take upper penta*/
    1138                 if (penta->IsOnSurface()) break;
    1139 
    1140                 /* get upper Penta*/
    1141                 penta=penta->GetUpperPenta();
    1142                 _assert_(penta->Id()!=this->id);
    1143         }
    1144 
    1145         /*return output*/
    1146         return penta;
    1147 }
    1148 /*}}}*/
    1149 Element*   Penta::GetUpperElement(void){/*{{{*/
    1150 
    1151         /*Output*/
    1152         Element* upper_element=this->GetUpperPenta();
    1153         return upper_element;
    1154 }
    1155 /*}}}*/
    11561129Penta*     Penta::GetUpperPenta(void){/*{{{*/
    11571130
    … …  
    13201293}
    13211294/*}}}*/
    1322 void       Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
    1323 
    1324         int  step,i;
     1295void       Penta::InputDepthAverageAtBase(int original_enum,int average_enum){/*{{{*/
     1296
     1297        IssmDouble  Jdet,value;
    13251298        IssmDouble  xyz_list[NUMVERTICES][3];
    1326         IssmDouble  Helem_list[NUMVERTICES];
    1327         IssmDouble  zeros_list[NUMVERTICES]={0.0};
    1328         IssmDouble  p0top1_list[NUMVERTICES];
    1329         Penta* penta=NULL;
    1330         Input* original_input=NULL;
    1331         Input* element_integrated_input=NULL;
    1332         Input* total_integrated_input=NULL;
    1333         Input* element_thickness_input=NULL;
    1334         Input* total_thickness_input=NULL;
    1335         Input* depth_averaged_input=NULL;
    1336 
    1337         /*recover parameters: */
     1299        IssmDouble  xyz_list_line[2][3];
     1300        IssmDouble  total[NUMVERTICES]       = {0.};
     1301        IssmDouble  intz[NUMVERTICES]        = {0.};
     1302        Input      *original_input           = NULL;
     1303        Input      *depth_averaged_input     = NULL;
    13381304
    13391305        /*Are we on the base? If not, return*/
    13401306        if(!IsOnBase()) return;
    13411307
    1342         /*OK, we are on bed. Initialize global inputs as 0*/
    1343         total_thickness_input =new PentaInput(ThicknessEnum,zeros_list,P1Enum);
    1344 
    13451308        /*Now follow all the upper element from the base to the surface to integrate the input*/
    1346         penta=this;
    1347         step =0;
     1309        Penta* penta = this;
     1310        int    step  = 0;
    13481311        for(;;){
    13491312
    1350                 /*Step1: Get original input (to be depth avegaged): */
    1351                 original_input=(Input*)penta->inputs->GetInput(enum_type);
    1352                 if(!original_input) _error_("could not find input with enum " << EnumToStringx(enum_type));
    1353 
    1354                 /*If first time, initialize total_integrated_input*/
    1355                 if (step==0){
    1356                         if (original_input->ObjectEnum()==PentaInputEnum)
    1357                          total_integrated_input=new PentaInput(average_enum_type,zeros_list,P1Enum);
    1358                         else if (original_input->ObjectEnum()==ControlInputEnum)
    1359                          total_integrated_input=new PentaInput(average_enum_type,zeros_list,P1Enum);
    1360                         else if (original_input->ObjectEnum()==DoubleInputEnum)
    1361                          total_integrated_input=new DoubleInput(average_enum_type,0.0);
    1362                         else{
    1363                          _error_("object " << EnumToStringx(original_input->ObjectEnum()) << " not supported yet");
    1364                         }
    1365                 }
     1313                /*Step1: Get original input (to be depth-avegaged): */
     1314                original_input=(Input*)penta->inputs->GetInput(original_enum);
     1315                if(!original_input) _error_("could not find input with enum " << EnumToStringx(original_enum));
    13661316
    13671317                /*Step2: Create element thickness input*/
    13681318                ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
    1369                 for(i=0;i<3;i++){
    1370                         Helem_list[i]=xyz_list[i+3][2]-xyz_list[i][2];
    1371                         Helem_list[i+3]=Helem_list[i];
    1372                 }
    1373                 element_thickness_input=new PentaInput(ThicknessEnum,Helem_list,P1Enum);
    1374 
    1375                 /*Step3: Vertically integrate A COPY of the original*/
    1376                 if(original_input->ObjectEnum()==PentaInputEnum){
    1377                         if(((PentaInput*)original_input)->interpolation_type==P0Enum){
    1378                                 original_input->GetInputValue(&p0top1_list[i]);
    1379                                 element_integrated_input= new  PentaInput(original_input->InstanceEnum(),p0top1_list,P1Enum);
    1380                                 element_integrated_input->VerticallyIntegrate(element_thickness_input);
     1319                for(int iv=0;iv<3;iv++){
     1320                        /*Get segment length*/
     1321                        for(int i=0;i<3;i++){
     1322                                xyz_list_line[0][i]=xyz_list[iv][i];
     1323                                xyz_list_line[1][i]=xyz_list[iv+3][i];
    13811324                        }
    1382                         else{
    1383                                 element_integrated_input= (Input*)original_input->copy();
    1384                                 element_integrated_input->VerticallyIntegrate(element_thickness_input);
     1325                        /*Integrate over edge*/
     1326                        Gauss* gauss=penta->NewGaussLine(iv,iv+3,3);
     1327                        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1328                                gauss->GaussPoint(ig);
     1329                                penta->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss);
     1330                                original_input->GetInputValue(&value,gauss);
     1331                                total[iv] += value*Jdet*gauss->weight;
     1332                                intz[iv]  += Jdet*gauss->weight;
    13851333                        }
    1386                 }
    1387                 else{
    1388                         element_integrated_input= (Input*)original_input->copy();
    1389                         element_integrated_input->VerticallyIntegrate(element_thickness_input);
    1390                 }
    1391 
    1392                 /*Add contributions to global inputs*/
    1393                 total_integrated_input->AXPY(element_integrated_input,1.0);
    1394                 total_thickness_input ->AXPY(element_thickness_input,1.0);
    1395 
    1396                 /*Clean up*/
    1397                 delete element_thickness_input;
    1398                 delete element_integrated_input;
     1334                        delete gauss;
     1335                }
    13991336
    14001337                /*Stop if we have reached the surface, else, take upper penta*/
    1401                 if (penta->IsOnSurface()) break;
     1338                if(penta->IsOnSurface()) break;
    14021339
    14031340                /* get upper Penta*/
    1404                 penta=penta->GetUpperPenta();
    1405                 _assert_(penta->Id()!=this->id);
    1406 
    1407                 /*increase couter*/
     1341                penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
    14081342                step++;
    14091343        }
    14101344
    1411         /*OK, now we only need to divide the depth integrated input by the total thickness!*/
    1412         depth_averaged_input=total_integrated_input->PointwiseDivide(total_thickness_input);
    1413         depth_averaged_input->ChangeEnum(average_enum_type);
    1414 
    1415         /*Clean up*/
    1416         delete total_thickness_input;
    1417         delete total_integrated_input;
     1345        /*Now we only need to divide the depth integrated input by the total thickness!*/
     1346        for(int iv=0;iv<3;iv++){
     1347                total[iv  ] = total[iv]/intz[iv];
     1348                total[iv+3] = total[iv];
     1349        }
     1350        switch(original_input->ObjectEnum()){
     1351                case PentaInputEnum:
     1352                case ControlInputEnum:
     1353                        depth_averaged_input=new PentaInput(average_enum,&total[0],P1Enum);
     1354                        break;
     1355                default:
     1356                        _error_("Interpolation " << EnumToStringx(original_input->ObjectEnum()) << " not supported yet");
     1357        }
    14181358
    14191359        /*Finally, add to inputs*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r21808 r21872  
    8585                Penta*         GetLowerPenta(void);
    8686                void           GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
    87                 Penta*         GetSurfacePenta(void);
    88                 Element*       GetUpperElement(void);
    8987                Penta*         GetUpperPenta(void);
    9088                void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r21808 r21872  
    7373                int         GetNumberOfVertices(void);
    7474                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
    75                 Element*    GetUpperElement(void){_error_("not implemented yet");};
    7675                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid){_error_("not implemented yet");};
    7776                void        GetVerticesCoordinates(IssmDouble** pxyz_list);

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r21808 r21872  
    7979                int         GetNumberOfVertices(void);
    8080                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
    81                 Element*    GetUpperElement(void){_error_("not implemented yet");};
    8281                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid){_error_("not implemented yet");};
    8382                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r21808 r21872  
    8484                int         GetNumberOfVertices(void);
    8585                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
    86                 Element*    GetUpperElement(void){_error_("not implemented yet");};
    8786                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
    8887                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);

  • issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp

    r21521 r21872  
    157157}
    158158/*}}}*/
    159 void BoolInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    160         /*square of a bool is the bool itself: */
    161         *psquaremin=value;
    162 }
    163 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h

    r21748 r21872  
    6969                void Scale(IssmDouble scale_factor);
    7070                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    71                 void SquareMin(IssmDouble* psquaremin, Parameters* parameters);
    72                 void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");};
    7371                /*}}}*/
    7472

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp

    r21416 r21872  
    282282        this->values->AXPY(gradient,scalar);
    283283}/*}}}*/
    284 void ControlInput::VerticallyIntegrate(Input* thickness_input){/*{{{*/
    285         values->VerticallyIntegrate(thickness_input);
    286 }/*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h

    r21748 r21872  
    8585                void SetGradient(Input* gradient_in);
    8686                void SetInput(Input* in_input);
    87                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
    8887                void UpdateValue(IssmDouble scalar);
    89                 void VerticallyIntegrate(Input* thickness_input);
    9088                /*}}}*/
    9189

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h

    r21748 r21872  
    8080                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    8181                void SetGradient(Input* gradient_in){_error_("not implemented yet");};
    82                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
    8382                void UpdateValue(IssmDouble scalar){_error_("not implemented yet");};
    84                 void VerticallyIntegrate(Input* thickness_input){_error_("not implemented yet");};
    8583                /*}}}*/
    8684

  • issm/trunk-jpl/src/c/classes/Inputs/DoubleArrayInput.h

    r21748 r21872  
    7171                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    7272                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    73                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
    74                 void VerticallyIntegrate(Input* thickness_input){_error_("not implemented yet");};
    7573                /*}}}*/
    7674

  • issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp

    r20827 r21872  
    245245}
    246246/*}}}*/
    247 void DoubleInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    248 
    249         /*square min of a IssmDouble is the square of the IssmDouble itself: */
    250         *psquaremin=pow(value,2);
    251 }
    252 /*}}}*/
    253 void DoubleInput::VerticallyIntegrate(Input* thickness_input){/*{{{*/
    254 
    255         /*Intermediaries*/
    256         IssmDouble thickness_value;
    257 
    258         /*Check that input provided is a thickness*/
    259         if (thickness_input->InstanceEnum()!=ThicknessEnum) _error_("Input provided is not a Thickness (enum_type is " << EnumToStringx(thickness_input->InstanceEnum()) << ")");
    260 
    261         /*vertically integrate depending on type:*/
    262         switch(thickness_input->ObjectEnum()){
    263 
    264                 case PentaInputEnum:
    265                         thickness_input->GetInputAverage(&thickness_value);
    266                         this->value=this->value*thickness_value;
    267                         return;
    268 
    269                 default:
    270                         _error_("not implemented yet");
    271         }
    272 }
    273 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h

    r21748 r21872  
    7272                void Scale(IssmDouble scale_factor);
    7373                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    74                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    75                 void VerticallyIntegrate(Input* thickness_input);
    7674                /*}}}*/
    7775

  • issm/trunk-jpl/src/c/classes/Inputs/Input.h

    r21748 r21872  
    5353                virtual void   Scale(IssmDouble scale_factor)=0;
    5454                virtual void   Set(IssmDouble setvalue)=0;
    55                 virtual void   SquareMin(IssmDouble* psquaremin,Parameters* parameters)=0;
    56                 virtual void   VerticallyIntegrate(Input* thickness_input)=0;
    5755
    5856                virtual int  GetResultArraySize(void)=0;

  • issm/trunk-jpl/src/c/classes/Inputs/IntInput.cpp

    r21521 r21872  
    161161}
    162162/*}}}*/
    163 void IntInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    164 
    165         /*square min of an integer is the square of the integer itself: */
    166         *psquaremin=pow((IssmDouble)value,2);
    167 }
    168 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/IntInput.h

    r21748 r21872  
    7373                void Scale(IssmDouble scale_factor);
    7474                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    75                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    76                 void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");};
    7775                /*}}}*/
    7876

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp

    r21586 r21872  
    407407}
    408408/*}}}*/
    409 void PentaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    410 
    411         int        numnodes=this->NumberofNodes(this->interpolation_type);
    412         IssmDouble squaremin;
    413 
    414         /*Now, figure out minimum of valuescopy: */
    415         squaremin=pow(this->values[0],2);
    416         for(int i=1;i<numnodes;i++){
    417                 if(pow(this->values[i],2)<squaremin)squaremin=pow(this->values[i],2);
    418         }
    419         /*Assign output pointers:*/
    420         *psquaremin=squaremin;
    421 }
    422 /*}}}*/
    423 void PentaInput::VerticallyIntegrate(Input* thickness_input){/*{{{*/
    424 
    425         IssmDouble thickness;
    426         IssmDouble value=0.;
    427 
    428         /*Check that input provided is a thickness*/
    429         if (thickness_input->InstanceEnum()!=ThicknessEnum) _error_("Input provided is not a Thickness (enum_type is " << EnumToStringx(thickness_input->InstanceEnum()) << ")");
    430 
    431         /*vertically integrate depending on type (and use P1 interpolation from now on)*/
    432         switch(this->interpolation_type){
    433                 case P1Enum:
    434                 case P1bubbleEnum:
    435                 case P1xP2Enum:
    436                 case P1xP3Enum:
    437                 case P2Enum:
    438                           {
    439                                 this->interpolation_type=P1Enum;
    440                                 GaussPenta *gauss=new GaussPenta();
    441                                 for(int iv=0;iv<3;iv++){
    442                                         gauss->GaussVertex(iv);
    443                                         thickness_input->GetInputValue(&thickness,gauss);
    444                                         this->values[iv]=0.5*(this->values[iv]+this->values[iv+3]) * thickness;
    445                                         this->values[iv+3]=this->values[iv];
    446                                 }
    447                                 delete gauss;
    448                                 return;
    449                           }
    450                 default:
    451                         _error_("not supported yet for type "<<EnumToStringx(this->interpolation_type));
    452         }
    453 }
    454 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h

    r21748 r21872  
    7373                Input* SpawnTriaInput(int index1,int index2,int index3);
    7474                Input* SpawnSegInput(int index1,int index2);
    75                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    76                 void VerticallyIntegrate(Input* thickness_input);
    7775
    7876};

  • issm/trunk-jpl/src/c/classes/Inputs/SegInput.h

    r21748 r21872  
    7474                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    7575                void Set(IssmDouble setvalue){_error_("not implemented yet");};
    76                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");};
    77                 void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");};
    7876
    7977};

  • issm/trunk-jpl/src/c/classes/Inputs/TetraInput.cpp

    r20827 r21872  
    418418}
    419419/*}}}*/
    420 void TetraInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    421 
    422         int        numnodes=this->NumberofNodes(this->interpolation_type);
    423         IssmDouble squaremin;
    424 
    425         /*Now, figure out minimum of valuescopy: */
    426         squaremin=pow(this->values[0],2);
    427         for(int i=1;i<numnodes;i++){
    428                 if(pow(this->values[i],2)<squaremin)squaremin=pow(this->values[i],2);
    429         }
    430         /*Assign output pointers:*/
    431         *psquaremin=squaremin;
    432 }
    433 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h

    r21748 r21872  
    7474                void Scale(IssmDouble scale_factor);
    7575                void Set(IssmDouble setvalue);
    76                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    77                 void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");};
    7876
    7977};

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

    r20827 r21872  
    535535}
    536536/*}}}*/
    537 void TransientInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    538 
    539         IssmDouble time;
    540 
    541         /*First, recover current time from parameters: */
    542         parameters->FindParam(&time,TimeEnum);
    543 
    544    /*Retrieve interpolated values for this time step: */
    545         Input* input=GetTimeInput(time);
    546 
    547         /*Call input function*/
    548         input->SquareMin(psquaremin,parameters);
    549 
    550         delete input;
    551 
    552 }
    553 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h

    r21748 r21872  
    8080                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    8181                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    82                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    83                 void VerticallyIntegrate(Input* thickness_forcing){_error_("not supported yet");};
    8482                /*}}}*/
    8583

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r20827 r21872  
    451451}
    452452/*}}}*/
    453 void TriaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/
    454 
    455         int        numnodes=this->NumberofNodes(this->interpolation_type);
    456         IssmDouble squaremin;
    457 
    458         /*Now, figure out minimum of valuescopy: */
    459         squaremin=pow(this->values[0],2);
    460         for(int i=1;i<numnodes;i++){
    461                 if(pow(this->values[i],2)<squaremin)squaremin=pow(this->values[i],2);
    462         }
    463         /*Assign output pointers:*/
    464         *psquaremin=squaremin;
    465 }
    466 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r21748 r21872  
    7474                void Scale(IssmDouble scale_factor);
    7575                void Set(IssmDouble setvalue);
    76                 void SquareMin(IssmDouble* psquaremin,Parameters* parameters);
    77                 void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");};
    7876
    7977};
Note: See TracChangeset for help on using the changeset viewer.