Changeset 17361


Ignore:
Timestamp:
02/27/14 10:15:38 (11 years ago)
Author:
bdef
Message:

CHG:transferring epl thickness computation to the analysis

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

Legend:

Unmodified
Added
Removed

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

    r17356 r17361  
    410410                break;
    411411        case 1:
    412 
    413                 element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum);
     412                /* _assert_(input) */
     413                /* get input */
     414
     415                        element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum);
    414416                element->GetInputValue(&sed_head,gauss,SedimentHeadEnum);
    415417                element->GetInputValue(&epl_head,gauss,EplHeadEnum);
    … …  
    481483        return transfer;
    482484}/*}}}*/
     485
     486void HydrologyDCEfficientAnalysis::ComputeEPLThickness(FemModel* femmodel){/*{{{*/
     487
     488        bool        active_element;
     489        int         meshtype;
     490        IssmDouble  dt,A,B;
     491        IssmDouble  EPLgrad2;
     492        IssmDouble  EPL_N;
     493
     494        femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
     495
     496        for(int j=0;j<femmodel->elements->Size();j++){
     497               
     498                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
     499               
     500                switch(meshtype){
     501                case Mesh2DhorizontalEnum:
     502                        if(!element->IsOnBed()) return;                 
     503                        B = element->GetMaterialParameter(MaterialsRheologyBbarEnum);
     504                        break;
     505                case Mesh3DEnum:
     506                        B = element->GetMaterialParameter(MaterialsRheologyBEnum);
     507                        break;
     508                default:
     509                _error_("not Implemented Yet");
     510                }
     511                       
     512                int         numnodes = element->GetNumberOfNodes();
     513                IssmDouble  thickness[numnodes];
     514                IssmDouble  eplhead[numnodes];
     515                IssmDouble  epl_slopeX[numnodes],epl_slopeY[numnodes];
     516                IssmDouble  old_thickness[numnodes];
     517                IssmDouble  ice_thickness[numnodes],bed[numnodes];
     518
     519                Input*  active_element_input=element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);               
     520                active_element_input->GetInputValue(&active_element);
     521                element->FindParam(&dt,TimesteppingTimeStepEnum);
     522       
     523                /*For now, assuming just one way to compute EPL thickness*/
     524                IssmDouble gravity          = element->GetMaterialParameter(ConstantsGEnum);
     525                IssmDouble rho_water        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     526                IssmDouble rho_ice          = element->GetMaterialParameter(MaterialsRhoIceEnum);
     527                IssmDouble n                =   element->GetMaterialParameter(MaterialsRheologyNEnum);
     528                IssmDouble latentheat       = element->GetMaterialParameter(MaterialsLatentheatEnum);
     529                IssmDouble epl_conductivity = element->GetMaterialParameter(HydrologydcEplConductivityEnum);
     530                IssmDouble init_thick       =   element->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);
     531               
     532                A=pow(B,-n);
     533               
     534                element->GetInputListOnVertices(&eplhead[0],EplHeadEnum);
     535                element->GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);
     536                element->GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
     537                element->GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
     538                element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
     539                element->GetInputListOnVertices(&bed[0],BedEnum);
     540                       
     541                if(!active_element){
     542                       
     543                        /*Keeping thickness to initial value if EPL is not active*/
     544                        for(int i=0;i<numnodes;i++){
     545                                thickness[i]=init_thick;
     546                        }
     547                }
     548                else{
     549                        for(int i=0;i<numnodes;i++){
     550                               
     551                                /*Compute first the effective pressure in the EPL*/
     552                                EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
     553                                if(EPL_N<0.0)EPL_N=0.0;
     554                                /*Get then the square of the gradient of EPL heads*/
     555                                EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i]+epl_slopeY[i]*epl_slopeY[i]);
     556                               
     557                                /*And proceed to the real thing*/
     558                                thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
     559                               
     560                                /*Take care of otherthikening*/
     561                                if(thickness[i]>10.0*init_thick){
     562                                        thickness[i] = 10.0*init_thick;
     563                                }
     564                        }
     565                }
     566                element->AddInput(HydrologydcEplThicknessEnum,thickness,P1Enum);
     567        }
     568}
     569/*}}}*/
    483570void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    484571        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.

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

    r17349 r17361  
    3737                IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss);
    3838                IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss);
     39                void ComputeEPLThickness(FemModel* femmodel);
    3940};
    4041#endif

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

    r17353 r17361  
    286286                virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;
    287287                virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
    288                 virtual void ComputeEPLThickness(void)=0;
     288                //              virtual void ComputeEPLThickness(void)=0;
    289289
    290290                virtual void   MigrateGroundingLine(IssmDouble* sheet_ungrounding)=0;

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

    r17355 r17361  
    42204220}
    42214221/*}}}*/
    4222 /*FUNCTION Penta::ComputeEPLThickness{{{*/
    4223 void  Penta::ComputeEPLThickness(void){
    4224 
    4225         int         i;
    4226         const int   numdof   = NDOF1 *NUMVERTICES;
    4227         const int   numdof2d = NDOF1 *NUMVERTICES2D;
    4228         bool        isefficientlayer;
    4229         IssmDouble  n,A,dt,init_thick;
    4230         IssmDouble  rho_water,rho_ice;
    4231         IssmDouble  gravity,latentheat,EPLgrad;
    4232         IssmDouble  EPL_N,epl_conductivity;
    4233         IssmDouble  activeEpl[numdof],thickness[numdof];
    4234         IssmDouble  eplhead[numdof], old_thickness[numdof];
    4235         IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof];
    4236         IssmDouble  ice_thickness[numdof],bed[numdof];
    4237         Penta       *penta = NULL;
    4238         /*If not on bed, return*/
    4239         if (!IsOnBed())return;
    4240 
    4241         /*Get the flag to know if the efficient layer is present*/
    4242         this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    4243         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4244 
    4245         if(isefficientlayer){
    4246                 /*For now, assuming just one way to compute EPL thickness*/
    4247                 rho_water        = matpar->GetRhoWater();
    4248                 rho_ice          = matpar->GetRhoIce();
    4249                 gravity          = matpar->GetG();
    4250                 latentheat       = matpar->GetLatentHeat();
    4251                 epl_conductivity = matpar->GetEplConductivity();
    4252                 init_thick       = matpar->GetEplInitialThickness();
    4253                 n                = material->GetN();
    4254                 A                = material->GetA();
     4222/* /\*FUNCTION Penta::ComputeEPLThickness{{{*\/ */
     4223/* void  Penta::ComputeEPLThickness(void){ */
     4224
     4225/*      int         i; */
     4226/*      const int   numdof   = NDOF1 *NUMVERTICES; */
     4227/*      const int   numdof2d = NDOF1 *NUMVERTICES2D; */
     4228/*      bool        isefficientlayer; */
     4229/*      IssmDouble  n,A,dt,init_thick; */
     4230/*      IssmDouble  rho_water,rho_ice; */
     4231/*      IssmDouble  gravity,latentheat,EPLgrad; */
     4232/*      IssmDouble  EPL_N,epl_conductivity; */
     4233/*      IssmDouble  activeEpl[numdof],thickness[numdof]; */
     4234/*      IssmDouble  eplhead[numdof], old_thickness[numdof]; */
     4235/*      IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof]; */
     4236/*      IssmDouble  ice_thickness[numdof],bed[numdof]; */
     4237/*      Penta       *penta = NULL; */
     4238/*      /\*If not on bed, return*\/ */
     4239/*      if (!IsOnBed())return; */
     4240
     4241/*      /\*Get the flag to know if the efficient layer is present*\/ */
     4242/*      this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */
     4243/*      this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */
     4244
     4245/*      if(isefficientlayer){ */
     4246/*              /\*For now, assuming just one way to compute EPL thickness*\/ */
     4247/*              rho_water        = matpar->GetRhoWater(); */
     4248/*              rho_ice          = matpar->GetRhoIce(); */
     4249/*              gravity          = matpar->GetG(); */
     4250/*              latentheat       = matpar->GetLatentHeat(); */
     4251/*              epl_conductivity = matpar->GetEplConductivity(); */
     4252/*              init_thick       = matpar->GetEplInitialThickness(); */
     4253/*              n                = material->GetN(); */
     4254/*              A                = material->GetA(); */
    42554255               
    4256                 GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum);
    4257                 GetInputListOnVertices(&eplhead[0],EplHeadEnum);
    4258                 GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);
    4259                 GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
    4260                 GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
    4261                 GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
    4262                 GetInputListOnVertices(&bed[0],BedEnum);
     4256/*              GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); */
     4257/*              GetInputListOnVertices(&eplhead[0],EplHeadEnum); */
     4258/*              GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);  */
     4259/*              GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */
     4260/*              GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */
     4261/*              GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */
     4262/*              GetInputListOnVertices(&bed[0],BedEnum); */
    42634263               
    4264                 for(int i=0;i<numdof2d;i++){
    4265                         /*Keeping thickness to 1 if EPL is not active*/
    4266                         if(activeEpl[i]==0.0){
    4267                                 thickness[i]=init_thick;
    4268                                 thickness[i+numdof2d]=thickness[i];
    4269                         }
    4270                         else{
    4271 
    4272                                 /*Compute first the effective pressure in the EPL*/
    4273                                 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
    4274                                 if(EPL_N<0.0)EPL_N=0.0;
    4275                                 /*Get then the gradient of EPL heads*/
    4276                                 EPLgrad = epl_slopeX[i]+epl_slopeY[i];
     4264/*              for(int i=0;i<numdof2d;i++){ */
     4265/*                      /\*Keeping thickness to 1 if EPL is not active*\/ */
     4266/*                      if(activeEpl[i]==0.0){ */
     4267/*                              thickness[i]=init_thick; */
     4268/*                              thickness[i+numdof2d]=thickness[i]; */
     4269/*                      } */
     4270/*                      else{ */
     4271
     4272/*                              /\*Compute first the effective pressure in the EPL*\/ */
     4273/*                              EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); */
     4274/*                              if(EPL_N<0.0)EPL_N=0.0; */
     4275/*                              /\*Get then the gradient of EPL heads*\/ */
     4276/*                              EPLgrad = epl_slopeX[i]+epl_slopeY[i]; */
    42774277                               
    4278                                 /*And proceed to the real thing*/
    4279                                 thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
    4280                                 thickness[i+numdof2d]=thickness[i];
    4281                         }
    4282                 }
    4283                 penta=this;
    4284                 for(;;){
    4285 
    4286                         /*Add input to the element: */                 
    4287                         penta->inputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));
     4278/*                              /\*And proceed to the real thing*\/ */
     4279/*                              thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); */
     4280/*                              thickness[i+numdof2d]=thickness[i]; */
     4281/*                      } */
     4282/*              } */
     4283/*              penta=this; */
     4284/*              for(;;){ */
     4285
     4286/*                      /\*Add input to the element: *\/                         */
     4287/*                      penta->inputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */
    42884288                       
    4289                         /*Stop if we have reached the surface*/
    4290                         if (penta->IsOnSurface()) break;
     4289/*                      /\*Stop if we have reached the surface*\/ */
     4290/*                      if (penta->IsOnSurface()) break; */
    42914291                       
    4292                         /* get upper Penta*/
    4293                         penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
    4294                 }
    4295         }
    4296 }
    4297 /*}}}*/
     4292/*                      /\* get upper Penta*\/ */
     4293/*                      penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); */
     4294/*              } */
     4295/*      } */
     4296/* } */
     4297/* /\*}}}*\/ */
    42984298
    42994299/*FUNCTION Penta::MigrateGroundingLine{{{*/

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

    r17355 r17361  
    244244                void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    245245                void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    246                 void           ComputeEPLThickness(void);
     246                //              void           ComputeEPLThickness(void);
    247247
    248248                void           UpdateConstraintsExtrudeFromBase(void);

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

    r17353 r17361  
    153153                void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
    154154                void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
    155                 void    ComputeEPLThickness(void){_error_("not implemented yet");};
     155                //              void    ComputeEPLThickness(void){_error_("not implemented yet");};
     156               
    156157                void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type){_error_("not implemented yet");};
    157158                void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum){_error_("not implemented yet");};

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

    r17355 r17361  
    45464546}
    45474547/*}}}*/
    4548 /*FUNCTION Tria::ComputeEPLThickness{{{*/
    4549 void  Tria::ComputeEPLThickness(void){
    4550 
    4551         const int   numdof         = NDOF1 *NUMVERTICES;
    4552         bool        isefficientlayer;
    4553         bool        active_element;
    4554         IssmDouble  n,A,dt,init_thick;
    4555         IssmDouble  rho_water,rho_ice;
    4556         IssmDouble  gravity,latentheat,EPLgrad2;
    4557         IssmDouble  EPL_N,epl_conductivity;
    4558         IssmDouble  thickness[numdof];
    4559         IssmDouble  eplhead[numdof];
    4560         IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof];
    4561         IssmDouble  old_thickness[numdof];
    4562         IssmDouble  ice_thickness[numdof],bed[numdof];
    4563 
    4564         Input* active_element_input=NULL;
    4565 
    4566         /*Get the flag to know if the efficient layer is present*/
    4567         this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    4568         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4569 
    4570         if(isefficientlayer){
    4571 
    4572                 /*For now, assuming just one way to compute EPL thickness*/
    4573                 rho_water        = matpar->GetRhoWater();
    4574                 rho_ice          = matpar->GetRhoIce();
    4575                 gravity          = matpar->GetG();
    4576                 latentheat       = matpar->GetLatentHeat();
    4577                 epl_conductivity = matpar->GetEplConductivity();
    4578                 init_thick       = matpar->GetEplInitialThickness();
    4579                 n                = material->GetN();
    4580                 A                = material->GetAbar();
     4548//* *FUNCTION Tria::ComputeEPLThickness{{{*\/ */
     4549/* void  Tria::ComputeEPLThickness(void){ */
     4550
     4551/*      const int   numdof         = NDOF1 *NUMVERTICES; */
     4552/*      bool        isefficientlayer; */
     4553/*      bool        active_element; */
     4554/*      IssmDouble  n,A,dt,init_thick; */
     4555/*      IssmDouble  rho_water,rho_ice; */
     4556/*      IssmDouble  gravity,latentheat,EPLgrad2; */
     4557/*      IssmDouble  EPL_N,epl_conductivity; */
     4558/*      IssmDouble  thickness[numdof]; */
     4559/*      IssmDouble  eplhead[numdof]; */
     4560/*      IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof]; */
     4561/*      IssmDouble  old_thickness[numdof]; */
     4562/*      IssmDouble  ice_thickness[numdof],bed[numdof]; */
     4563
     4564/*      Input* active_element_input=NULL; */
     4565
     4566/*      /\*Get the flag to know if the efficient layer is present*\/ */
     4567/*      this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */
     4568/*      this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */
     4569
     4570/*      if(isefficientlayer){ */
     4571
     4572/*              /\*For now, assuming just one way to compute EPL thickness*\/ */
     4573/*              rho_water        = matpar->GetRhoWater(); */
     4574/*              rho_ice          = matpar->GetRhoIce(); */
     4575/*              gravity          = matpar->GetG(); */
     4576/*              latentheat       = matpar->GetLatentHeat(); */
     4577/*              epl_conductivity = matpar->GetEplConductivity(); */
     4578/*              init_thick       = matpar->GetEplInitialThickness(); */
     4579/*              n                = material->GetN(); */
     4580/*              A                = material->GetAbar(); */
    45814581               
    4582                 active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
    4583                 active_element_input->GetInputValue(&active_element);
     4582/*              active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */
     4583/*              active_element_input->GetInputValue(&active_element); */
    45844584                       
    4585                 GetInputListOnVertices(&eplhead[0],EplHeadEnum);
    4586                 GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);
    4587                 GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
    4588                 GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
    4589                 GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
    4590                 GetInputListOnVertices(&bed[0],BedEnum);
     4585/*              GetInputListOnVertices(&eplhead[0],EplHeadEnum); */
     4586/*              GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);  */
     4587/*              GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */
     4588/*              GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */
     4589/*              GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */
     4590/*              GetInputListOnVertices(&bed[0],BedEnum); */
    45914591               
    4592                 if(!active_element){
    4593                         /*Keeping thickness to initial value if EPL is not active*/
    4594                         for(int i=0;i<numdof;i++){
    4595                                 thickness[i]=init_thick;
    4596                         }
    4597                 }
    4598                 else{
    4599                         for(int i=0;i<numdof;i++){
     4592/*              if(!active_element){ */
     4593/*                      /\*Keeping thickness to initial value if EPL is not active*\/ */
     4594/*                      for(int i=0;i<numdof;i++){ */
     4595/*                              thickness[i]=init_thick; */
     4596/*                      } */
     4597/*              } */
     4598/*              else{ */
     4599/*                      for(int i=0;i<numdof;i++){ */
    46004600                               
    4601                                 /*Compute first the effective pressure in the EPL*/
    4602                                 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
    4603                                 if(EPL_N<0.0)EPL_N=0.0;
    4604                                 /*Get then the square of the gradient of EPL heads*/
    4605                                 EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i]+epl_slopeY[i]*epl_slopeY[i]);
     4601/*                              /\*Compute first the effective pressure in the EPL*\/ */
     4602/*                              EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); */
     4603/*                              if(EPL_N<0.0)EPL_N=0.0; */
     4604/*                              /\*Get then the square of the gradient of EPL heads*\/ */
     4605/*                              EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i]+epl_slopeY[i]*epl_slopeY[i]); */
    46064606                               
    4607                                 /*And proceed to the real thing*/
    4608                                 thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
     4607/*                              /\*And proceed to the real thing*\/ */
     4608/*                              thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); */
    46094609                                       
    4610                                 /*Take care of otherthikening*/
    4611                                 if(thickness[i]>10.0*init_thick){
    4612                                         thickness[i] = 10.0*init_thick;
    4613                                 }
    4614                         }
    4615                 }
    4616                 this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));
    4617         }
    4618 }
     4610/*                              /\*Take care of otherthikening*\/ */
     4611/*                              if(thickness[i]>10.0*init_thick){ */
     4612/*                                      thickness[i] = 10.0*init_thick; */
     4613/*                              } */
     4614/*                      } */
     4615/*              } */
     4616/*              this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */
     4617/*      } */
     4618/* } */
    46194619/*}}}*/
    46204620

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

    r17355 r17361  
    253253
    254254                void           CreateHydrologyWaterVelocityInput(void);
     255               
    255256                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    256257                void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    257258                void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    258                 void           ComputeEPLThickness(void);
     259                //              void           ComputeEPLThickness(void);
    259260                /*}}}*/
    260261

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

    r17353 r17361  
    14071407}
    14081408/*}}}*/
    1409 void FemModel::HydrologyEPLThicknessx(void){ /*{{{*/
    1410 
    1411         for (int i=0;i<elements->Size();i++){
    1412                 Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    1413                 element->ComputeEPLThickness();
    1414         }
    1415 }
     1409/* void FemModel::HydrologyEPLThicknessx(void){ /\*{{{*\/ */
     1410
     1411/*      for (int i=0;i<elements->Size();i++){ */
     1412/*              Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); */
     1413/*              element->ComputeEPLThickness(); */
     1414/*      } */
     1415/* } */
    14161416/*}}}*/
    14171417void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/

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

    r17353 r17361  
    9696                void UpdateConstraintsExtrudeFromTopx();
    9797                void HydrologyEPLupdateDomainx(void);
    98                 void HydrologyEPLThicknessx(void);
     98                //              void HydrologyEPLThicknessx(void);
    9999                void UpdateConstraintsL2ProjectionEPLx(void);
    100100};

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp

    r17275 r17361  
    639639                if(zigzag_counter>penalty_lock){
    640640                        unstable=0;
    641                         active=1;
     641                        active=0;
    642642                }
    643643        }

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

    r17357 r17361  
    3333        Vector<IssmDouble>* df=NULL;
    3434
     35        HydrologyDCInefficientAnalysis* inefanalysis = NULL;
     36        HydrologyDCEfficientAnalysis* effanalysis = NULL;
     37       
    3538        bool       sedconverged,eplconverged,hydroconverged;
    3639        bool       isefficientlayer;
    … …  
    6063
    6164        if(isefficientlayer) {
     65                inefanalysis = new HydrologyDCInefficientAnalysis();
     66                effanalysis = new HydrologyDCEfficientAnalysis();
    6267                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
    6368                GetSolutionFromInputsx(&ug_epl,femmodel);
    64                
    6569                /*Initialize the element mask*/
    66                 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
    67                 analysis->ElementizeEplMask(femmodel);
    68                 delete analysis;
     70                inefanalysis->ElementizeEplMask(femmodel);
    6971        }
    7072        /*The real computation starts here, outermost loop is on the two layer system*/
    … …  
    115117                                        if(num_unstable_constraints==0) sedconverged = true;
    116118                                        if (sedcount>=hydro_maxiter){
     119                                                //sedconverged = true;
    117120                                                _error_("   maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded");
    118121                                        }
    … …  
    126129                        /*Update EPL mask*/
    127130                        if(isefficientlayer){
    128                                 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
    129                                 analysis->ElementizeEplMask(femmodel);
    130                                 delete analysis;
     131                                inefanalysis->ElementizeEplMask(femmodel);
    131132                        }
    132133                        sedconverged=false;
    … …  
    180181                                solutionsequence_linear(femmodel);
    181182                                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
    182                                 femmodel->HydrologyEPLThicknessx();
    183                                
     183
     184                                effanalysis->ComputeEPLThickness(femmodel);
     185
    184186                                //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
    185187                                femmodel->HydrologyEPLupdateDomainx();
    186                                
    187                                 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
    188                                 analysis->ElementizeEplMask(femmodel);
    189                                 delete analysis;
     188                                inefanalysis->ElementizeEplMask(femmodel);
    190189                                /* }}} */
    191190                                       
    … …  
    281280        delete uf_epl;
    282281        delete uf_epl_sub_iter;
    283        
     282        delete inefanalysis;
     283        delete effanalysis;
    284284}
Note: See TracChangeset for help on using the changeset viewer.