Changeset 19094


Ignore:
Timestamp:
02/06/15 16:44:07 (10 years ago)
Author:
bdef
Message:

CHG: Adding node capbility for inputupdatefromvector and changes in the zigzag counter in epl

Location:
issm/trunk-jpl/src/c
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

    r19092 r19094  
    672672        int         eplflip_lock;
    673673        int         numnodes      =basalelement->GetNumberOfNodes();
    674         IssmDouble  new_active;
    675         IssmDouble  recurent;
    676674        IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
    677675        IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
    … …  
    700698        for(i=0;i<numnodes;i++){
    701699
    702                 /*Dealing with the initial value to define zigzaging, preceding iteration the first time we see the node and current evolving vector after*/
    703                 recurence->GetValue(&recurent,basalelement->nodes[i]->Sid());
    704                 if (recurent>0)vec_mask->GetValue(&old_active[i],basalelement->nodes[i]->Sid());
    705                 new_active=old_active[i];
    706                 recurence->SetValue(basalelement->nodes[i]->Sid(),1.,ADD_VAL);
     700                /*If node is now closed bring its thickness back to initial*/
     701                if (old_active[i]==0.){
     702                        epl_thickness[i]=init_thick;
     703                }
    707704
    708705                /*Now starting to look at the activations*/
    709706                if(residual[i]>0.){
    710707                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    711                         new_active=1.0;
     708                        if(old_active[i]==0.)   recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    712709                }
    713710                /*If mask was already one, keep one or colapse*/
    714711                else if(old_active[i]>0.){
    715712                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    716                         new_active=1.0;
    717713                        /*If epl thickness gets under colapse thickness, close the layer*/
    718714                        if(epl_thickness[i]<colapse_thick){
    719715                                vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
    720                                 new_active=0.0;
     716                                recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    721717                        }
    722718                }
    … …  
    728724                                if(sedhead[j] == sedheadmin){
    729725                                        vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
     726                                        if(old_active[i]==0.)   recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
    730727                                }
    731728                        }
    732                 }
    733                 /*Now checking for nodes zigzaging between open and close*/
    734                 if (old_active[i] != new_active){
    735                         eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
    736                         /*If node changed too much of state, fix it to it's last known state*/
    737                         if(eplzigzag_counter[basalelement->nodes[i]->Lid()]>eplflip_lock & eplflip_lock!=0){
    738                                 vec_mask->SetValue(basalelement->nodes[i]->Sid(),old_active[i],INS_VAL);
    739                                 new_active=old_active[i];
    740                         }
    741                 }
    742                         /*If node is now closed bring its thickness back to initial*/
    743                 if (new_active==0.){
    744                         epl_thickness[i]=init_thick;
    745729                }
    746730        }

  • issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

    r19033 r19094  
    556556
    557557        Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
    558         GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
     558        GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
    559559       
    560560        /* set NaN on elements intersected by zero levelset */

  • issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.cpp

    r18930 r19094  
    417417
    418418        Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
    419         GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
     419        GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
    420420       
    421421        /* set distance on elements intersected by zero levelset */

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

    r19019 r19094  
    862862}
    863863/*}}}*/
    864 void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/*{{{*/
     864/* void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/\*{{{*\/ */
     865
     866/*      /\*Fetch number vertices for this element and allocate arrays*\/ */
     867/*      int numvertices = this->GetNumberOfVertices(); */
     868/*      int*        vertexpidlist = xNew<int>(numvertices); */
     869/*      IssmDouble* values        = xNew<IssmDouble>(numvertices); */
     870
     871/*      /\*Fill in values*\/ */
     872/*      this->GetVertexPidList(vertexpidlist); */
     873/*      this->GetInputListOnVertices(values,input_enum); */
     874/*      vector->SetValues(numvertices,vertexpidlist,values,INS_VAL); */
     875
     876/*      /\*Clean up*\/ */
     877/*      xDelete<int>(vertexpidlist); */
     878/*      xDelete<IssmDouble>(values); */
     879
     880/* } */
     881/* /\*}}}*\/ */
     882void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type){/*{{{*/
    865883
    866884        /*Fetch number vertices for this element and allocate arrays*/
    867         int numvertices = this->GetNumberOfVertices();
    868         int*        vertexpidlist = xNew<int>(numvertices);
    869         IssmDouble* values        = xNew<IssmDouble>(numvertices);
    870 
    871         /*Fill in values*/
    872         this->GetVertexPidList(vertexpidlist);
    873         this->GetInputListOnVertices(values,input_enum);
    874         vector->SetValues(numvertices,vertexpidlist,values,INS_VAL);
    875 
     885        int         numvertices = this->GetNumberOfVertices();
     886        int         numnodes    = this->GetNumberOfNodes();
     887        int*        doflist     = NULL;
     888        IssmDouble* values      = NULL;
     889
     890        switch(type){
     891        case VertexPIdEnum:
     892                doflist = xNew<int>(numvertices);
     893                values = xNew<IssmDouble>(numvertices);
     894                /*Fill in values*/
     895                this->GetVertexPidList(doflist);
     896                this->GetInputListOnVertices(values,input_enum);
     897                vector->SetValues(numvertices,doflist,values,INS_VAL);
     898                break;
     899        case VertexSIdEnum:
     900                doflist = xNew<int>(numvertices);
     901                values = xNew<IssmDouble>(numvertices);
     902                /*Fill in values*/
     903                this->GetVerticesSidList(doflist);
     904                this->GetInputListOnVertices(values,input_enum);
     905                vector->SetValues(numvertices,doflist,values,INS_VAL);
     906                break;
     907        case NodesEnum:
     908                doflist = xNew<int>(numnodes);
     909                values = xNew<IssmDouble>(numnodes);
     910                /*Fill in values*/
     911                this->GetInputListOnNodes(values,input_enum);
     912                this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     913                vector->SetValues(numnodes,doflist,values,INS_VAL);
     914                break;
     915        case NodeSIdEnum:
     916                doflist = xNew<int>(numnodes);
     917                values = xNew<IssmDouble>(numnodes);
     918                /*Fill in values*/
     919                this->GetNodesSidList(doflist);
     920                this->GetInputListOnNodes(values,input_enum);
     921                vector->SetValues(numnodes,doflist,values,INS_VAL);
     922                break;
     923        default:
     924                _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
     925        }
     926       
    876927        /*Clean up*/
    877         xDelete<int>(vertexpidlist);
     928        xDelete<int>(doflist);
    878929        xDelete<IssmDouble>(values);
    879930

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

    r19061 r19094  
    9696                void               GetNodesSidList(int* sidlist);
    9797                void               GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
    98                 void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
     98                void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum, int type);
    9999                void                 GetVertexPidList(int* pidlist);
    100100                void               GetVerticesConnectivityList(int* connectivitylist);

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

    r19085 r19094  
    19431943
    19441944        case NodesEnum:
    1945 
    19461945                /*Get number of nodes and dof list: */
    19471946                numnodes = this->NumberofNodes(this->element_type);
    … …  
    19571956
    19581957        case NodeSIdEnum:
    1959 
    19601958                /*Get number of nodes and dof list: */
    19611959                numnodes = this->NumberofNodes(this->element_type);

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

    r19091 r19094  
    611611
    612612        /*get vertex vectors for bed and thickness: */
    613         GetVectorFromInputsx(&surface  ,this, SurfaceEnum,VertexEnum);
    614         GetVectorFromInputsx(&bed      ,this, BaseEnum,    VertexEnum);
     613        GetVectorFromInputsx(&surface  ,this, SurfaceEnum,VertexPIdEnum);
     614        GetVectorFromInputsx(&bed      ,this, BaseEnum,   VertexPIdEnum);
    615615
    616616        /*Allocate vector*/
    … …  
    18311831
    18321832                        /*this response was scaled. pick up the response from the inputs: */
    1833                         GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexEnum);
     1833                        GetVectorFromInputsx(&vertex_response,this, StringToEnumx(root),VertexPIdEnum);
    18341834
    18351835                        /*Now, average it onto the partition nodes: */
    … …  
    19141914        Vector<IssmDouble>* mask                                                        = NULL;
    19151915        Vector<IssmDouble>* recurence                           = NULL;
     1916        Vector<IssmDouble>* active                                              = NULL;
    19161917        IssmDouble*         serial_mask                         = NULL;
    1917         Vector<IssmDouble>* active                                              = NULL;
     1918        IssmDouble*         serial_rec                          = NULL;
    19181919        IssmDouble*         serial_active                       = NULL;
     1920        IssmDouble*         old_active        = NULL;
    19191921        int*                eplzigzag_counter = NULL;
    1920                
     1922        int                 eplflip_lock;
     1923       
    19211924        HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
    19221925        HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
     1926
    19231927        /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
    19241928        mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
    19251929        recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
    19261930        this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
     1931        this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
     1932        GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
     1933       
    19271934        for (int i=0;i<elements->Size();i++){
    19281935                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    19291936                effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element);
    19301937        }
     1938        /*check for changes and increment zigzag counter, change the mask if necessary*/
     1939        recurence->Assemble();
     1940        serial_rec=recurence->ToMPISerial();
     1941        for (int i=0;i<nodes->Size();i++){
     1942                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
     1943                if(serial_rec[node->Sid()]==1.)eplzigzag_counter[node->Lid()] ++;
     1944                if(eplzigzag_counter[node->Lid()]>eplflip_lock & eplflip_lock!=0){
     1945                        mask->SetValue(node->Sid(),old_active[node->Sid()],INS_VAL);
     1946                }
     1947        }
     1948
     1949       
    19311950        this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
    19321951        /*Assemble and serialize*/
    19331952        mask->Assemble();
    1934         serial_mask=mask->ToMPISerial();
     1953        serial_mask=mask->ToMPISerial();       
     1954       
    19351955        xDelete<int>(eplzigzag_counter);
    19361956        delete mask;

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

    r18521 r19094  
    77#include "../../toolkits/toolkits.h"
    88
    9 void GetVectorFromInputsx( Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type){
     9void GetVectorFromInputsx(Vector<IssmDouble>** pvector,FemModel* femmodel,int name,int type){
    1010
    1111        int i;
    1212        Vector<IssmDouble>* vector=NULL;
    1313
    14         if(type==VertexEnum){
    15 
    16                 /*Allocate vector*/
     14        switch(type){
     15        case VertexPIdEnum: case VertexSIdEnum:
    1716                vector=new Vector<IssmDouble>(femmodel->vertices->NumberOfVertices());
    18 
    19                 /*Look up in elements*/
    20                 for(i=0;i<femmodel->elements->Size();i++){
    21                         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    22                         element->GetVectorFromInputs(vector,name);
    23                 }
     17                break;
     18        case NodesEnum:case NodeSIdEnum:
     19                vector=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
     20                break;
     21        default:
     22                        _error_("vector type: " << EnumToStringx(type) << " not supported yet!");
    2423        }
    25         else{
    26                 _error_("vector type: " << EnumToStringx(type) << " not supported yet!");
     24        /*Look up in elements*/
     25        for(i=0;i<femmodel->elements->Size();i++){
     26                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     27                element->GetVectorFromInputs(vector,name,type);
    2728        }
    2829

  • issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp

    r18909 r19094  
    142142        for(i=0;i<elements->Size();i++){
    143143                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    144                 element->GetVectorFromInputs(vec_phi,MaskGroundediceLevelsetEnum);
     144                element->GetVectorFromInputs(vec_phi,MaskGroundediceLevelsetEnum,VertexPIdEnum);
    145145        }
    146146        vec_phi->Assemble();

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp

    r18654 r19094  
    4444        Vector<IssmDouble>* vel_old       = NULL;
    4545        Vector<IssmDouble>* vel_old_local = NULL;
    46         GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
    47         GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
     46        GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum);
     47        GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum);
    4848
    4949        while(true){
    … …  
    7070                        /*Update solution*/
    7171                        InputUpdateFromSolutionx(femmodel,ug); delete ug;
    72                         GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexEnum);
     72                        GetVectorFromInputsx(&vel,femmodel,VxEnum,VertexPIdEnum);
    7373                        /*Check for convergence*/
    7474                        Vector<IssmDouble>* dvel_local=vel_old_local->Duplicate(); vel_old_local->Copy(dvel_local); dvel_local->AYPX(vel,-1.0);
    … …  
    101101                /*Update solution*/
    102102                InputUpdateFromSolutionx(femmodel,pug); delete pug;
    103                 GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexEnum);
     103                GetVectorFromInputsx(&pug,femmodel,PressureEnum,VertexPIdEnum);
    104104
    105105                /*Check for convergence*/

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp

    r17678 r19094  
    4141        Vector<IssmDouble>* vx     = NULL;
    4242        Vector<IssmDouble>* vx_old = NULL;
    43         GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
     43        GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum);
    4444
    4545        while(true){
    … …  
    6969                analysis->InputUpdateFromSolutionFSXTH_d(  femmodel->elements,femmodel->parameters);
    7070                analysis->InputUpdateFromSolutionFSXTH_tau(femmodel->elements,femmodel->parameters);
    71                 GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum);
     71                GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexPIdEnum);
    7272
    7373                /*Check for convergence*/
Note: See TracChangeset for help on using the changeset viewer.