Changeset 15249


Ignore:
Timestamp:
06/13/13 08:36:55 (12 years ago)
Author:
bdef
Message:

CHG: Changes in the way we deal with the EPL mask

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

Legend:

Unmodified
Added
Removed

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

    r15223 r15249  
    4141                virtual void   CreatePVector(Vector<IssmDouble>* pf)=0;
    4242                virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
    43                 virtual void   BasisIntegral(Vector<IssmDouble>* basisg)=0;
    4443                virtual void   GetSolutionFromInputs(Vector<IssmDouble>* solution)=0;
    4544                virtual int    GetNodeIndex(Node* node)=0;
    … …  
    136135                virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
    137136                virtual void UpdateConstraints(void)=0;
     137                virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;
     138                virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
    138139                #endif
    139140};

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

    r15223 r15249  
    802802        /*return output*/
    803803        return penta;
    804 }
    805 /*}}}*/
    806 /*FUNCTION Penta::BasisIntegral {{{*/
    807 void Penta::BasisIntegral(Vector<IssmDouble>* basisg){
    808 
    809         /*Constants*/
    810         const int    numdof=NDOF1*NUMVERTICES;
    811 
    812         /*Intermediaries */
    813         IssmDouble Jdet;
    814         IssmDouble xyz_list[NUMVERTICES][3];
    815         IssmDouble basis[numdof];
    816         IssmDouble basisint[numdof]={0.};
    817         int       *doflist=NULL;
    818         GaussPenta* gauss=NULL;
    819 
    820         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    821         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    822 
    823         /* Start looping on the number of gaussian points: */
    824         gauss=new GaussPenta(2,2);
    825         for(int ig=gauss->begin();ig<gauss->end();ig++){
    826 
    827                 gauss->GaussPoint(ig);
    828 
    829                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    830                 GetNodalFunctionsP1(&basis[0], gauss);
    831 
    832                 for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];
    833         }
    834 
    835         basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);
    836 
    837         /*Clean up and return*/
    838         delete gauss;
    839         xDelete<int>(doflist);
    840804}
    841805/*}}}*/
    … …  
    21662130                                name==InversionVyObsEnum ||
    21672131                                name==InversionVzObsEnum ||
    2168                                 name==BasisIntegralEnum ||
    21692132                                name==TemperatureEnum ||
    21702133                                name==EnthalpyEnum ||
    … …  
    92989261}
    92999262/*}}}*/
     9263
     9264/*FUNCTION Penta::HydrologyEPLGetActive {{{*/
     9265void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
     9266        _error_("Hydrological stuff not suported in Penta");
     9267}
     9268/*}}}*/
     9269/*FUNCTION Penta::HydrologyEPLGetMask{{{*/
     9270void  Penta::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){
     9271        _error_("Hydrological stuff not suported in Penta");
     9272}
     9273/*}}}*/
    93009274#endif

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

    r15223 r15249  
    263263                void           InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong);
    264264                void           InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong);
    265                 void           BasisIntegral(Vector<IssmDouble>* gbasis);
    266265                void             GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solutiong);
    267266                void             GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solutiong);
    … …  
    307306                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    308307                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
     308                void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
     309                void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    309310                #endif
    310311                #ifdef _HAVE_THERMAL_

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

    r15234 r15249  
    2626/*FUNCTION Tria::Tria(){{{*/
    2727Tria::Tria(){
    28 
     28       
    2929        int i;
    3030
    … …  
    10991099        _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
    11001100        return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
    1101 }
    1102 /*}}}*/
    1103 /*FUNCTION Tria::BasisIntegral {{{*/
    1104 void Tria::BasisIntegral(Vector<IssmDouble>* basisg){
    1105 
    1106         /*Constants*/
    1107         const int    numdof=NDOF1*NUMVERTICES;
    1108 
    1109         /*Intermediaries */
    1110         IssmDouble Jdet;
    1111         IssmDouble xyz_list[NUMVERTICES][3];
    1112         IssmDouble basis[numdof];
    1113         IssmDouble basisint[numdof]={0.};
    1114         int       *doflist=NULL;
    1115         GaussTria* gauss=NULL;
    1116 
    1117         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    1118         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    1119 
    1120         /* Start looping on the number of gaussian points: */
    1121         gauss=new GaussTria(2);
    1122         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1123 
    1124                 gauss->GaussPoint(ig);
    1125 
    1126                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    1127                 GetNodalFunctions(&basis[0], gauss);
    1128 
    1129                 for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];
    1130         }
    1131 
    1132         basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);
    1133 
    1134         /*Clean up and return*/
    1135         delete gauss;
    1136         xDelete<int>(doflist);
    11371101}
    11381102/*}}}*/
    … …  
    18181782                        break;
    18191783                case HydrologyDCEfficientAnalysisEnum:
    1820                         InputUpdateFromSolutionHydrologyDCEfficient(solution);
     1784                        InputUpdateFromSolutionOneDof(solution,EplHeadEnum);
    18211785                        break;
    18221786                #endif
    … …  
    18491813void  Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){
    18501814
    1851         const int numdof          = NDOF1*NUMVERTICES;
    1852 
    1853         int*      doflist=NULL;
    1854         IssmDouble    values[numdof];
     1815        const int  numdof         = NDOF1*NUMVERTICES;
     1816        int*       doflist        = NULL;
     1817        IssmDouble values[numdof];
    18551818
    18561819        /*Get dof list: */
    … …  
    19811944                for(int i=0;i<numdof;i++){
    19821945                        values[i]=vector[doflist[i]];
     1946                        if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     1947                }
     1948                /*Add input to the element: */
     1949                this->inputs->AddInput(new TriaP1Input(name,values));
     1950
     1951                /*Free ressources:*/
     1952                xDelete<int>(doflist);
     1953                return;
     1954
     1955        case NodeSIdEnum:
     1956                for(int i=0;i<NUMVERTICES;i++){
     1957                        values[i]=vector[nodes[i]->Sid()];
    19831958                        if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    19841959                }
    … …  
    21392114                                name==HydrologydcMaskEplactiveEnum ||
    21402115                                name==WaterTransferEnum ||
    2141                                 name==BasisIntegralEnum ||
    21422116                                name==QmuVxEnum ||
    21432117                                name==QmuVyEnum ||
    … …  
    28282802        /*Skip if water element*/
    28292803        if(IsOnWater()) return;
    2830 
    2831         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    2832         switch(analysis_type){
    2833                 #ifdef _HAVE_HYDROLOGY_
    2834                 case HydrologyDCEfficientAnalysisEnum:
    2835                         this->HydrologyUpdateEplDomain();
    2836                         break;
    2837                 #endif
    2838         }
    28392804
    28402805}
    … …  
    58295794/*FUNCTION Tria::AllActive{{{*/
    58305795bool Tria::AllActive(void){
    5831 
     5796       
    58325797        /*Intermediaries*/
    58335798        const int  numnodes = NUMVERTICES;
    5834 
     5799       
    58355800        for(int i=0;i<numnodes;i++){
    58365801                if(!this->nodes[i]->IsActive()) return false;
    58375802        }
    5838 
     5803       
    58395804        return true;
    58405805}
    58415806/*}}}*/
     5807/*FUNCTION Tria::AnyActive{{{*/
     5808bool Tria::AnyActive(void){
     5809       
     5810        /*Intermediaries*/
     5811        const int  numnodes = NUMVERTICES;
     5812       
     5813        for(int i=0;i<numnodes;i++){
     5814                if(this->nodes[i]->IsActive()) return true;
     5815        }
     5816       
     5817        return false;
     5818}/*}}}*/
    58425819/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{*/
    58435820void Tria::CreateHydrologyWaterVelocityInput(void){
    … …  
    60766053        /*Check that all nodes are active, else return empty matrix*/
    60776054        if(!this->AllActive()) return NULL;
    6078 
    60796055        /*Initialize Element matrix*/
    60806056        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    … …  
    65006476}
    65016477/*}}}*/
    6502 /*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCEfficient{{{*/
    6503 void  Tria::InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution){
    6504 
    6505         /*Intermediaries*/
    6506         const int   numdof         = NDOF1 *NUMVERTICES;
    6507         int        *doflist        = NULL;
    6508         IssmDouble  values[numdof];
    6509 
    6510         /*Get dof list: */
    6511         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    6512 
    6513         /*Use the dof list to index into the solution vector: */
    6514         for(int i=0;i<numdof;i++){
    6515                 values[i]=solution[doflist[i]];
    6516                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    6517         }
    6518 
    6519         /*Add input to the element: */
    6520         this->inputs->AddInput(new TriaP1Input(EplHeadEnum,values));
    6521 
    6522         /*Free ressources:*/
    6523         xDelete<int>(doflist);
    6524 }
    6525 /*}}}*/
    65266478/*FUNCTION Tria::GetHydrologyDCInefficientHmax{{{*/
    65276479void  Tria::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){
    … …  
    65706522        IssmDouble sed_trans,sed_thick;
    65716523        IssmDouble leakage,h_max;
     6524        IssmDouble wh_trans;
    65726525        IssmDouble activeEpl[numdof];
    6573         IssmDouble wh_trans[numdof]={0.0};
    65746526        IssmDouble storing[numdof];
    65756527        IssmDouble epl_head[numdof],sed_head[numdof];
    65766528
    65776529        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    6578 
     6530       
    65796531        /*Get the flag to know if the efficient layer is present*/
    6580                 this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    6581 
     6532        this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     6533       
    65826534        if(isefficientlayer){
    65836535                /*Also get the flag to the transfer method*/
    … …  
    65926544
    65936545                        GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
    6594                         GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
     6546                        GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); 
    65956547                        GetInputListOnVertices(&epl_head[0],EplHeadEnum);
    65966548
    … …  
    66056557                                if(sed_head[i]>epl_head[i]){
    66066558                                        storing[i]=matpar->GetSedimentStoring();
    6607                                         wh_trans[i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
     6559                                        wh_trans=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
    66086560                                }
    66096561                                else{
    66106562                                        this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
    66116563                                        storing[i]=matpar->GetEplStoring();
    6612                                         wh_trans[i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
     6564                                        wh_trans=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
    66136565                                        if(sed_head[i]>h_max){
    6614                                                 wh_trans[i]=0.0;
     6566                                                wh_trans=0.0;
    66156567                                        }
    6616                                         if((sed_head[i]+wh_trans[i])>h_max){
    6617                                                 wh_trans[i]=h_max-sed_head[i];
     6568                                        if((sed_head[i]+wh_trans)>h_max){
     6569                                                wh_trans=h_max-sed_head[i];
    66186570                                        }
    66196571                                }
    6620                                 transfer->SetValue(doflist[i],wh_trans[i],INS_VAL);
     6572                                /*Assign output pointer*/
     6573                                transfer->SetValue(doflist[i],wh_trans,INS_VAL);
    66216574                        }
    66226575                        break;
    … …  
    66256578                }
    66266579        }
    6627         /*Assign output pointer*/
    6628         /*transfer->SetValues(numdof,doflist,&wh_trans[0],INS_VAL);*/
    6629 }
    6630 
    6631 /*}}}*/
    6632 /*FUNCTION Tria::HydrologyUpdateEplDomain{{{*/
    6633 void  Tria::HydrologyUpdateEplDomain(void){
     6580}
     6581/*}}}*/
     6582/*FUNCTION Tria::HydrologyEPLGetActive {{{*/
     6583void Tria::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
     6584
     6585        /*Constants*/
     6586        const int  numnodes         =NUMVERTICES;
     6587        int        flag;
     6588        IssmDouble active[numnodes];
     6589       
     6590        GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveEnum);
     6591               
     6592        for(int i=0;i<numnodes;i++) flag+=active[i];
     6593       
     6594        if(flag>0.){
     6595                for(int i=0;i<numnodes;i++) nodes[i]->Activate();
     6596        }
     6597       
     6598        else if(!this->AnyActive()){
     6599                for(int i=0;i<numnodes;i++) nodes[i]->Deactivate();
     6600               
     6601        }
     6602        else{
     6603                /*Do not do anything: at least one node is active for this element but this element is not solved for*/
     6604        }
     6605       
     6606}
     6607/*}}}*/
     6608/*FUNCTION Tria::HydrologyEPLGetMask{{{*/
     6609void  Tria::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){
    66346610
    66356611        /*Intermediaries*/
    6636         const int  numdof = NDOF1 *NUMVERTICES;
    6637         bool       isefficientlayer;
    6638         IssmDouble activeEpl[numdof];
    6639 
    6640         /*Get parameter and Inout list*/
    6641         this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    6642         GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
    6643 
    6644         if(isefficientlayer){
    6645                 for(int i=0;i<numdof;i++){
    6646                         if(activeEpl[i]==1.0)
    6647                          nodes[i]->Activate();
    6648                         else if (activeEpl[i]==0.0)
    6649                          nodes[i]->Deactivate();
    6650                         else
    6651                          _error_("activation value "<< activeEpl[i] <<" not supported");
     6612        int         i,j;
     6613        const int   numdof         = NDOF1 *NUMVERTICES;
     6614        IssmDouble  h_max;
     6615        IssmDouble  sedheadmin;
     6616        IssmDouble  old_active[numdof];
     6617        IssmDouble  sedhead[numdof];
     6618        IssmDouble  eplhead[numdof];
     6619
     6620
     6621        GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveEnum);   
     6622        GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
     6623        GetInputListOnVertices(&eplhead[0],EplHeadEnum);
     6624
     6625        sedheadmin=sedhead[0];
     6626        for(i=1;i<numdof;i++){
     6627                if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
     6628        }
     6629       
     6630        for(i=0;i<numdof;i++){
     6631                /*If mask was alread one, keep one*/
     6632                if(old_active[i]>0.){
     6633                        vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
     6634                }
     6635               
     6636                this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
     6637                /*Increase the size of the efficient system if needed*/
     6638                /*Increase is needed if the epl head reach the maximum value (sediment value for now)*/
     6639                if(eplhead[i]>=h_max){
     6640                        vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
     6641                        for(j=0;j<numdof;j++){
     6642
     6643                                /*Increase of the domain is on the downstream node in term of sediment head*/
     6644                                if(sedhead[j] == sedheadmin){
     6645                                        vec_mask->SetValue(nodes[j]->Sid(),1.,INS_VAL);
     6646                                        break;
     6647                                }
     6648                        }
    66526649                }
    66536650        }

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

    r15230 r15249  
    8989                int    Sid();
    9090                bool   IsOnBed();
    91                 bool   IsFloating(); 
     91                bool   IsFloating();
    9292                bool   IsNodeOnShelf();
    9393                bool   IsNodeOnShelfFromFlags(IssmDouble* flags);
    … …  
    196196                ElementVector* CreatePVectorSlope(void);
    197197                IssmDouble     GetArea(void);
    198                 void           BasisIntegral(Vector<IssmDouble>* gbasis);
    199198                int            GetElementType(void);
    200199                void             GetDofList(int** pdoflist,int approximation_enum,int setenum);
    … …  
    250249                void      GetSolutionFromInputsHydrologyDCEfficient(Vector<IssmDouble>* solution);
    251250                void    CreateHydrologyWaterVelocityInput(void);
    252                 void    HydrologyUpdateEplDomain(void);
    253251                void    UpdateConstraints(void);
    254252                void      InputUpdateFromSolutionHydrology(IssmDouble* solution);
    … …  
    259257                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    260258                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
     259                void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
     260                void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    261261                bool    AllActive(void);
     262                bool    AnyActive(void);
    262263                #endif
    263264                #ifdef _HAVE_BALANCED_

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

    r15237 r15249  
    15811581#endif
    15821582
    1583 void FemModel::BasisIntegralsx(void){ /*{{{*/
    1584 
    1585         Vector<IssmDouble>* basisg=NULL;
    1586 
    1587         /*Vector allocation*/
    1588         basisg=new Vector<IssmDouble>(nodes->NumberOfNodes());
    1589 
     1583void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/
     1584
     1585        Vector<IssmDouble> *mask          = NULL;
     1586        IssmDouble         *serial_mask   = NULL;
     1587        Vector<IssmDouble> *active        = NULL;
     1588        IssmDouble         *serial_active = NULL;
     1589
     1590        /*Step 1: update maks, the mask might be extended by residual and/or using downstream sediment head*/
     1591        mask=new Vector<IssmDouble>(nodes->NumberOfNodes());
    15901592        for (int i=0;i<elements->Size();i++){
    15911593                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    1592                 element->BasisIntegral(basisg);
    1593         }
    1594         /*Assemble*/
    1595         basisg->Assemble();
    1596 
    1597         /*Update Inputs*/
    1598         InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,basisg,BasisIntegralEnum,NodesEnum);
     1594                element->HydrologyEPLGetMask(mask);
     1595        }
     1596
     1597        /*Assemble and serialize*/
     1598        mask->Assemble();
     1599        serial_mask=mask->ToMPISerial();
     1600        delete mask;
     1601
     1602        /*Update Mask*/
     1603        InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,serial_mask,HydrologydcMaskEplactiveEnum,NodeSIdEnum);
     1604        xDelete<IssmDouble>(serial_mask);
     1605
     1606        /*Step 2: update node activity. If one element is connected to mask=1, all nodes are active*/
     1607        active=new Vector<IssmDouble>(nodes->NumberOfNodes());
     1608        for (int i=0;i<elements->Size();i++){
     1609                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1610                element->HydrologyEPLGetActive(active);
     1611        }
     1612
     1613        /*Assemble and serialize*/
     1614        active->Assemble();
     1615        serial_active=active->ToMPISerial();
     1616        delete active;
     1617
     1618        /*Update node activation accordingly*/
     1619        for (int i=0;i<nodes->Size();i++){
     1620                Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
     1621                if(serial_active[node->Sid()]==1.){
     1622                        node->Activate();
     1623                }
     1624                else{
     1625                        node->Deactivate();
     1626                }
     1627        }
     1628        xDelete<IssmDouble>(serial_active);
    15991629
    16001630}

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

    r15130 r15249  
    9999                void UpdateConstraintsx(void);
    100100                int  UpdateVertexPositionsx(void);
    101                 void BasisIntegralsx(void);             
     101                void ParEplMask(void);         
    102102                void HydrologyTransferx(void);
     103                void HydrologyEPLupdateDomainx(void);
    103104};
    104105

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

    r15234 r15249  
    4141        femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum);
    4242        femmodel->parameters->FindParam(&time,TimeEnum);
    43         femmodel->BasisIntegralsx();
    4443        hydro_maxiter=100;
    4544        hydrocount=1;
    … …  
    123122                                delete uf;
    124123                                Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
    125                                 delete old_uf; old_uf=uf->Duplicate();                 
     124                                delete old_uf; old_uf=uf->Duplicate();
    126125                                if(eplcount>1) delete ug_epl;
    127126                                delete Kff;delete pf;
    … …  
    129128                                Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
    130129                                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
     130                                femmodel->HydrologyEPLupdateDomainx();
    131131                                ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    132132                               
    … …  
    144144                                        InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum);
    145145                                        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
     146                                        femmodel->HydrologyEPLupdateDomainx();
    146147                                        break;
    147148                                }
Note: See TracChangeset for help on using the changeset viewer.