source: issm/oecreview/Archive/24684-25833/ISSM-25488-25489.diff

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 14.7 KB

  • ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp

     
    3535        /*Assign output pointers:*/
    3636        *pvector=vector;
    3737} /*}}}*/
     38void GetVectoronBaseFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type){ /*{{{*/
     39
     40        int i, domaintype;
     41        Vector<IssmDouble>* vector=NULL;
     42
     43        switch(type){
     44                case ElementSIdEnum:
     45                        vector=new Vector<IssmDouble>(femmodel->elements->NumberOfElements());
     46                        break;
     47                case VertexPIdEnum: case VertexSIdEnum:
     48                        vector=new Vector<IssmDouble>(femmodel->vertices->NumberOfVertices());
     49                        break;
     50                case NodesEnum:case NodeSIdEnum:
     51                        vector=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
     52                        break;
     53                default:
     54                        _error_("vector type: " << EnumToStringx(type) << " not supported yet!");
     55        }
     56
     57        /*Look up in elements*/
     58        for(i=0;i<femmodel->elements->Size();i++){
     59                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     60                element->FindParam(&domaintype,DomainTypeEnum);
     61                switch(domaintype){
     62                        case Domain2DhorizontalEnum:
     63                                element->GetVectorFromInputs(vector,name,type);
     64                                break;
     65                        case Domain3DEnum:
     66                                if(!element->IsOnBase()) continue;
     67                                element->GetVectorFromInputs(vector,name,type);
     68                                break;
     69                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     70                }
     71        }
     72        vector->Assemble();
     73        /*Assign output pointers:*/
     74        *pvector=vector;
     75} /*}}}*/
    3876void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time){/*{{{*/
    3977
    4078        int i;
    … …  
    78116        /*Assign output pointers:*/
    79117        *pvector=vector;
    80118}/*}}}*/
     119void GetVectoronBaseFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name, int type){/*{{{*/
     120
     121        /*output: */
     122        IssmDouble* vector=NULL;
     123
     124        /*intermediary: */
     125        Vector<IssmDouble>* vec_vector=NULL;
     126
     127        GetVectoronBaseFromInputsx(&vec_vector,femmodel,name,type);
     128        vector=vec_vector->ToMPISerial();
     129
     130        /*Free ressources:*/
     131        delete vec_vector;
     132
     133        /*Assign output pointers:*/
     134        *pvector=vector;
     135}/*}}}*/
    81136void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name){ /*{{{*/
    82137
    83138        int interpolation_type;
    84         /*this one is special: we don't specify the type, but let the nature of the inputs dictace. 
     139        /*this one is special: we don't specify the type, but let the nature of the inputs dictace.
    85140         * P0 -> ElementSIdEnum, P1 ->VertexSIdEnum: */
    86141
    87142        /*We go find the input of the first element, and query its interpolation type: */
    88143        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(0));
    89         Input* input=element->GetInput(name); 
     144        Input* input=element->GetInput(name);
    90145        if (!input) _error_("could not find input: " << name);
    91146
    92         interpolation_type=input->GetInputInterpolationType(); 
     147        interpolation_type=input->GetInputInterpolationType();
    93148        if(interpolation_type==P0Enum){
    94149                *pvector_size=femmodel->elements->NumberOfElements();
    95150                GetVectorFromInputsx(pvector,femmodel,name, ElementSIdEnum);

  • ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h

     
    11/*!\file:  GetVectorFromInputsx.h
    2  */ 
     2 */
    33
    44#ifndef _GETVECTORFROMINPUTSXX_H
    55#define _GETVECTORFROMINPUTSXX_H
    … …  
    77#include "../../classes/classes.h"
    88
    99/* local prototypes: */
    10 void GetVectorFromInputsx( Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type);
    11 void GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type);
    12 void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name);
    13 void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time);
     10void    GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type);
     11void    GetVectoronBaseFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type);
     12void  GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type,IssmDouble time);
     13void    GetVectorFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name,int type);
     14void    GetVectoronBaseFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name,int type);
     15void  GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name);
    1416
     17
    1518#endif  /* _GETVECTORFROMINPUTSXX_H */

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

     
    115115                                        /*Proceed now to heads computations*/
    116116                                        solutionsequence_hydro_nonlinear(femmodel);
    117117               /*If we have a sub-timestep we store the substep inputs in a transient input here*/
    118                                         femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
     118                                        femmodel->StackTransientInputonBasex(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
    119119                                }
    120120                                /*averaging the stack*/
    121                                 femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
     121                                femmodel->AverageTransientInputonBasex(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
    122122                                /*And reseting to global time*/
    123123                                femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
    124124                                femmodel->parameters->SetParam(global_time,TimeEnum);

  • ../trunk-jpl/src/c/classes/FemModel.h

     
    178178                void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
    179179                void InitTransientInputx(int* transientinput_enum,int numoutputs);
    180180                void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
     181                void StackTransientInputonBasex(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
    181182                void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method);
     183                void AverageTransientInputonBasex(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method);
    182184                void UpdateConstraintsx(void);
    183185                int  UpdateVertexPositionsx(void);
    184186

  • ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

     
    687687
    688688        bool     element_active;
    689689        Element* element=NULL;
    690 
    691         for(int i=0;i<femmodel->elements->Size();i++){
     690        int      elementssize=femmodel->elements->Size();
     691        for(int i=0;i<elementssize;i++){
    692692                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    693693
    694694                Input* input=element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(input);
    … …  
    747747
    748748        bool     element_active;
    749749        Element* element=NULL;
    750 
    751         for(int i=0;i<femmodel->elements->Size();i++){
     750        int elementssize = femmodel->elements->Size();
     751        for(int i=0;i<elementssize;i++){
    752752                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    753753
    754754                Input* input=element->GetInput(HydrologydcMaskThawedNodeEnum); _assert_(input);

  • ../trunk-jpl/src/c/classes/FemModel.cpp

     
    4848#include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
    4949#include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
    5050#include "../modules/NodalValuex/NodalValuex.h"
    51 #include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
    5251#include "../modules/AverageOntoPartitionx/AverageOntoPartitionx.h"
    5352/*}}}*/
    5453
    … …  
    45424541        int          npart;
    45434542
    45444543        /*retrieve partition vectors for responses that are scaled:*/
    4545         this->parameters->FindParam(&response_partitions,&response_partitions_num,NULL,NULL,QmuResponsePartitionsEnum); 
    4546         this->parameters->FindParam(&response_partitions_npart,NULL,NULL,QmuResponsePartitionsNpartEnum); 
    4547        
     4544        this->parameters->FindParam(&response_partitions,&response_partitions_num,NULL,NULL,QmuResponsePartitionsEnum);
     4545        this->parameters->FindParam(&response_partitions_npart,NULL,NULL,QmuResponsePartitionsNpartEnum);
     4546
    45484547        /*retrieve my_rank: */
    45494548        my_rank=IssmComm::GetRank();
    45504549
    … …  
    50205019        recurence=new Vector<IssmDouble>(numnodes);
    50215020        this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
    50225021        this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
    5023         GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
     5022        GetVectoronBaseFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
    50245023
    50255024        for (int i=0;i<elements->Size();i++){
    50265025                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    … …  
    50675066
    50685067        /*Update node activation accordingly*/
    50695068        int         counter  = 0; //this is probably not acurate but we are only interested in positivity
     5069        Element*    basalelement=NULL;
    50705070        for(int i=0;i<elements->Size();i++){
    50715071                Element    *element  = xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    5072                 int         numnodes = element->GetNumberOfNodes();
     5072                switch(domaintype){
     5073                        case Domain2DhorizontalEnum:
     5074                                basalelement = element;
     5075                                break;
     5076                        case Domain3DEnum:
     5077                                if(!element->IsOnBase()) continue;
     5078                                basalelement = element->SpawnBasalElement();
     5079                                break;
     5080                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     5081                }
     5082
     5083                int         numnodes = basalelement->GetNumberOfNodes();
    50735084                IssmDouble *base     = xNew<IssmDouble>(numnodes);
    5074                 element->GetInputListOnNodes(&base[0],BaseEnum);
     5085                basalelement->GetInputListOnNodes(&base[0],BaseEnum);
    50755086                for(int in=0;in<numnodes;in++){
    5076                         Node* node=element->GetNode(in);
     5087                        Node* node=basalelement->GetNode(in);
    50775088                        if(serial_active[node->Sid()]==1.){
    50785089                                node->Activate();
    50795090                                if(!node->IsClone()) counter++;
    … …  
    50845095                        }
    50855096                }
    50865097                xDelete<IssmDouble>(base);
     5098                if(domaintype!=Domain2DhorizontalEnum){
     5099                        basalelement->DeleteMaterials();
     5100                        delete basalelement;
     5101                }
    50875102        }
    50885103        xDelete<IssmDouble>(serial_active);
    50895104        delete effanalysis;
    … …  
    51335148        }
    51345149        /*for other cases we just grab the mask from the initialisation value*/
    51355150        else{
    5136                 GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
     5151                GetVectoronBaseFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
    51375152        }
    51385153        /*Update Mask and elementize*/
    51395154        InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
    … …  
    52805295        }
    52815296}
    52825297/*}}}*/
     5298void FemModel::StackTransientInputonBasex(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/
     5299
     5300        Element*   basalelement=NULL;
     5301        int      domaintype;
     5302        this->parameters->FindParam(&domaintype,DomainTypeEnum);
     5303
     5304        for(int i=0;i<numoutputs;i++){
     5305                if(input_enum[i]<0){
     5306                        _error_("Can't deal with non enum fields for result Stack");
     5307                }
     5308                else{
     5309                        for(int j=0;j<elements->Size();j++){
     5310                                /*Get the right transient input*/
     5311                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     5312                                /*Get basal element*/
     5313                                switch(domaintype){
     5314                                        case Domain2DhorizontalEnum:
     5315                                                break;
     5316                                        case Domain3DEnum:
     5317                                                if(!element->IsOnBase()) continue;
     5318                                                break;
     5319                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     5320                                }
     5321                                /*Get values and lid list*/
     5322                                TransientInput* transientinput = this->inputs->GetTransientInput(transientinput_enum[i]);
     5323                                const int   numvertices = element->GetNumberOfVertices();
     5324                                IssmDouble* values      = xNew<IssmDouble>(numvertices);
     5325                                int        *vertexlids  = xNew<int>(numvertices);
     5326                                switch(domaintype){
     5327                                        case Domain2DhorizontalEnum:
     5328                                                element->GetInputListOnVertices(&values[0],input_enum[i]); //this is the enum to stack
     5329                                                element->GetVerticesLidList(vertexlids);
     5330                                                transientinput->AddTriaTimeInput(subtime,numvertices,vertexlids,values,P1Enum);
     5331                                                break;
     5332                                        case Domain3DEnum:{
     5333                                                element->GetInputListOnVertices(&values[0],input_enum[i]); //this is the enum to stack
     5334                                                element->GetVerticesLidList(vertexlids);
     5335                                                Penta* penta=xDynamicCast<Penta*>(elements->GetObjectByOffset(j));
     5336                                                for(;;){
     5337                                                        transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum);
     5338                                                        if (penta->IsOnSurface()) break;
     5339                                                        penta=penta->GetUpperPenta(); _assert_(penta->Id()!=element->id);
     5340                                                        penta->GetVerticesLidList(vertexlids);
     5341                                                }
     5342                                                break;
     5343                                        }
     5344                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     5345                                }
     5346                                xDelete<IssmDouble>(values);
     5347                                xDelete<int>(vertexlids);
     5348                        }
     5349                }
     5350        }
     5351}
     5352/*}}}*/
    52835353void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs, int averaging_method){ /*{{{*/
    52845354
    52855355        for(int i=0;i<numoutputs;i++){
    … …  
    52895359                }
    52905360        }
    52915361}/*}}}*/
     5362void FemModel::AverageTransientInputonBasex(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs, int averaging_method){ /*{{{*/
     5363
     5364        int      domaintype;
     5365        this->parameters->FindParam(&domaintype,DomainTypeEnum);
     5366
     5367        for(int i=0;i<numoutputs;i++){
     5368                for(int j=0;j<this->elements->Size();j++){
     5369                        Element* element = xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     5370                        /*Get basal element*/
     5371                        switch(domaintype){
     5372                                case Domain2DhorizontalEnum:
     5373                                        element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);;
     5374                                        break;
     5375                                case Domain3DEnum:
     5376                                        if(!element->IsOnBase()) continue;
     5377                                        element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);
     5378                                        break;
     5379                                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     5380                        }
     5381                }
     5382        }
     5383}/*}}}*/
    52925384#ifdef _HAVE_JAVASCRIPT_
    52935385FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
    52945386        /*configuration: */
Note: See TracBrowser for help on using the repository browser.