Changeset 27649


Ignore:
Timestamp:
03/18/23 13:04:20 (2 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on making SHAKTI work in 3d, not there yet :(

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

Legend:

Unmodified
Added
Removed

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

    r27470 r27649  
    5353        int M,N;
    5454        int *segments = NULL;
    55         iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
     55        if(iomodel->domaintype==Domain3DEnum){
     56                iomodel->FetchData(&segments,&M,&N,"md.mesh.segments2d");
     57        }
     58        else if(iomodel->domaintype==Domain2DhorizontalEnum){
     59                iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
     60        }
     61        else{
     62                _error_("mesh type not supported yet");
     63        }
    5664
    5765        /*Check that the size seem right*/
    5866        _assert_(N==3); _assert_(M>=3);
     67
    5968        for(int i=0;i<M;i++){
    6069                if(iomodel->my_elements[segments[i*3+2]-1]){
    … …  
    169178ElementMatrix* HydrologyShaktiAnalysis::CreateKMatrix(Element* element){/*{{{*/
    170179
     180        /* Check if ice in element */
     181        if(element->IsAllFloating() || !element->IsIceInElement()) return NULL;
     182        if(!element->IsOnBase()) return NULL;
     183        Element* basalelement = element->SpawnBasalElement();
     184
    171185        /*Intermediaries */
    172186        IssmDouble Jdet;
    … …  
    175189
    176190        /*Fetch number of nodes and dof for this finite element*/
    177         int numnodes = element->GetNumberOfNodes();
     191        int numnodes = basalelement->GetNumberOfNodes();
    178192
    179193        /*Initialize Element vector and other vectors*/
    180         ElementMatrix* Ke     = element->NewElementMatrix();
     194        ElementMatrix* Ke     = basalelement->NewElementMatrix();
    181195        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    182196        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    183197
    184198        /*Retrieve all inputs and parameters*/
    185         element->GetVerticesCoordinates(&xyz_list);
     199        basalelement->GetVerticesCoordinates(&xyz_list);
    186200
    187201        /*Get conductivity from inputs*/
    188         IssmDouble conductivity = GetConductivity(element);
     202        IssmDouble conductivity = GetConductivity(basalelement);
    189203
    190204        /*Get englacial storage coefficient*/
    191205        IssmDouble storage,dt;
    192         element->FindParam(&storage,HydrologyStorageEnum);
    193         element->FindParam(&dt,TimesteppingTimeStepEnum);
    194 
    195         /*Get all inputs and parameters*/
    196         element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum);
    197         element->FindParam(&rho_ice,MaterialsRhoIceEnum);
    198         element->FindParam(&g,ConstantsGEnum);
    199         Input* B_input = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    200         Input* n_input = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
    201         Input* gap_input = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
    202         Input* thickness_input = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
    203         Input* head_input = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
    204         Input* base_input = element->GetInput(BaseEnum);                      _assert_(base_input);
     206        basalelement->FindParam(&storage,HydrologyStorageEnum);
     207        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     208
     209        /*Get all inputs and parameters*/
     210        basalelement->FindParam(&rho_water,MaterialsRhoFreshwaterEnum);
     211        basalelement->FindParam(&rho_ice,MaterialsRhoIceEnum);
     212        basalelement->FindParam(&g,ConstantsGEnum);
     213        Input* B_input = basalelement->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
     214        Input* n_input = basalelement->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     215        Input* gap_input = basalelement->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
     216        Input* thickness_input = basalelement->GetInput(ThicknessEnum);                  _assert_(thickness_input);
     217        Input* head_input = basalelement->GetInput(HydrologyHeadEnum);              _assert_(head_input);
     218        Input* base_input = basalelement->GetInput(BaseEnum);                      _assert_(base_input);
    205219
    206220        /* Start  looping on the number of gaussian points: */
    207         Gauss* gauss=element->NewGauss(1);
     221        Gauss* gauss=basalelement->NewGauss(1);
    208222        while(gauss->next()){
    209223
    210                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    211                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    212                 element->NodalFunctions(basis,gauss);
    213 
    214                 base_input->GetInputValue(&bed,gauss);
    215                 thickness_input->GetInputValue(&thickness,gauss);
    216                 gap_input->GetInputValue(&gap,gauss);
    217                 head_input->GetInputValue(&head,gauss);
    218 
    219                 /*Get ice A parameter*/
    220                 B_input->GetInputValue(&B,gauss);
    221                 n_input->GetInputValue(&n,gauss);
    222                 A=pow(B,-n);
    223 
    224                 /*Get water and ice pressures*/
    225                 IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.);
    226                 IssmDouble pressure_water = rho_water*g*(head-bed);
    227                 if(pressure_water>pressure_ice) pressure_water = pressure_ice;
    228 
    229 
     224                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     225                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     226                basalelement->NodalFunctions(basis,gauss);
     227
     228                base_input->GetInputValue(&bed,gauss);
     229                thickness_input->GetInputValue(&thickness,gauss);
     230                gap_input->GetInputValue(&gap,gauss);
     231                head_input->GetInputValue(&head,gauss);
     232
     233                /*Get ice A parameter*/
     234                B_input->GetInputValue(&B,gauss);
     235                n_input->GetInputValue(&n,gauss);
     236                A=pow(B,-n);
     237
     238                /*Get water and ice pressures*/
     239                IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.);
     240                IssmDouble pressure_water = rho_water*g*(head-bed);
     241                if(pressure_water>pressure_ice) pressure_water = pressure_ice;
    230242
    231243                for(int i=0;i<numnodes;i++){
    … …  
    233245                                Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
    234246                                  + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]
    235                                         +gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j];
     247                                  +gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j];
    236248                        }
    237249                }
    … …  
    243255        xDelete<IssmDouble>(dbasis);
    244256        delete gauss;
     257        if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
    245258        return Ke;
    246259}/*}}}*/
    … …  
    248261
    249262        /*Skip if water or ice shelf element*/
    250         if(element->IsAllFloating()) return NULL;
     263        if(element->IsAllFloating() || !element->IsIceInElement()) return NULL;
     264        if(!element->IsOnBase()) return NULL;
     265        Element* basalelement = element->SpawnBasalElement();
    251266
    252267        /*Intermediaries */
    … …  
    259274
    260275        /*Fetch number of nodes and dof for this finite element*/
    261         int numnodes = element->GetNumberOfNodes();
     276        int numnodes = basalelement->GetNumberOfNodes();
    262277
    263278        /*Initialize Element vector and other vectors*/
    264         ElementVector* pe    = element->NewElementVector();
     279        ElementVector* pe    = basalelement->NewElementVector();
    265280        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    266281
    267282        /*Retrieve all inputs and parameters*/
    268         element->GetVerticesCoordinates(&xyz_list);
    269         IssmDouble  latentheat      = element->FindParam(MaterialsLatentheatEnum);
    270         IssmDouble  g               = element->FindParam(ConstantsGEnum);
    271         IssmDouble  rho_ice         = element->FindParam(MaterialsRhoIceEnum);
    272         IssmDouble  rho_water       = element->FindParam(MaterialsRhoFreshwaterEnum);
    273         Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
    274         Input* head_input           = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
    275         Input* gap_input            = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
    276         Input* thickness_input      = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
    277         Input* base_input           = element->GetInput(BaseEnum);                       _assert_(base_input);
    278         Input* B_input              = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    279         Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
    280         Input* englacial_input      = element->GetInput(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
    281         Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
    282         Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    283    Input* headold_input        = element->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
     283        basalelement->GetVerticesCoordinates(&xyz_list);
     284        IssmDouble  latentheat      = basalelement->FindParam(MaterialsLatentheatEnum);
     285        IssmDouble  g               = basalelement->FindParam(ConstantsGEnum);
     286        IssmDouble  rho_ice         = basalelement->FindParam(MaterialsRhoIceEnum);
     287        IssmDouble  rho_water       = basalelement->FindParam(MaterialsRhoFreshwaterEnum);
     288        Input* geothermalflux_input = basalelement->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
     289        Input* head_input           = basalelement->GetInput(HydrologyHeadEnum);              _assert_(head_input);
     290        Input* gap_input            = basalelement->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
     291        Input* thickness_input      = basalelement->GetInput(ThicknessEnum);                  _assert_(thickness_input);
     292        Input* base_input           = basalelement->GetInput(BaseEnum);                       _assert_(base_input);
     293        Input* B_input              = basalelement->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
     294        Input* n_input              = basalelement->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     295        Input* englacial_input      = basalelement->GetInput(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
     296        Input* lr_input             = basalelement->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
     297        Input* br_input             = basalelement->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
     298   Input* headold_input        = basalelement->GetInput(HydrologyHeadOldEnum);           _assert_(headold_input);
    284299
    285300        /*Get conductivity from inputs*/
    286         IssmDouble conductivity = GetConductivity(element);
     301        IssmDouble conductivity = GetConductivity(basalelement);
    287302
    288303        /*Get englacial storage coefficient*/
    289304        IssmDouble storage,dt;
    290    element->FindParam(&storage,HydrologyStorageEnum);
    291    element->FindParam(&dt,TimesteppingTimeStepEnum);
    292 
    293         /*Build friction element, needed later: */
    294         Friction* friction=new Friction(element,2);
     305   basalelement->FindParam(&storage,HydrologyStorageEnum);
     306   basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     307
     308        /*Build friction basalelement, needed later: */
     309        Friction* friction=new Friction(basalelement,2);
    295310
    296311        /* Start  looping on the number of gaussian points: */
    297         Gauss* gauss=element->NewGauss(2);
     312        Gauss* gauss=basalelement->NewGauss(2);
    298313        while(gauss->next()){
    299314
    300                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    301                 element->NodalFunctions(basis,gauss);
     315                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     316                basalelement->NodalFunctions(basis,gauss);
    302317                geothermalflux_input->GetInputValue(&G,gauss);
    303318                base_input->GetInputValue(&bed,gauss);
    … …  
    338353
    339354                /*Compute change in sensible heat due to changes in pressure melting point*/
    340                 dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
     355                dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
    341356                dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
    342357
    343         meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
    344 
    345                   for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    346                   (
    347                     meltrate*(1/rho_water-1/rho_ice)
    348                     +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
    349                     +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
    350                     -beta*sqrt(vx*vx+vy*vy)
    351                     +ieb
    352                     +storage*head_old/dt
    353                     )*basis[i];
    354 
    355        
    356         }
     358                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
     359
     360                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
     361                (
     362                  meltrate*(1/rho_water-1/rho_ice)
     363                  +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap
     364                  +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap
     365                  -beta*sqrt(vx*vx+vy*vy)
     366                  +ieb
     367                  +storage*head_old/dt
     368                )*basis[i];
     369
     370        }
     371
    357372        /*Clean up and return*/
    358373        xDelete<IssmDouble>(xyz_list);
    … …  
    360375        delete friction;
    361376        delete gauss;
     377        if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
    362378        return pe;
    363379}/*}}}*/
    … …  
    369385}/*}}}*/
    370386void           HydrologyShaktiAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
     387
     388        /*Only update if on base*/
     389        if(!element->IsOnBase()) return;
    371390
    372391        /*Intermediary*/
    … …  
    411430
    412431        /*Add input to the element: */
    413         element->AddInput(HydrologyHeadEnum,values,element->GetElementType());
     432        element->AddBasalInput(HydrologyHeadEnum,values,element->GetElementType());
    414433
    415434        /*Update reynolds number according to new solution*/
    … …  
    424443
    425444        IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU;
    426         element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum);
     445        element->AddBasalInput(HydrologyReynoldsEnum,&reynolds,P0Enum);
    427446
    428447   /*Compute new effective pressure*/
    … …  
    438457}/*}}}*/
    439458void           HydrologyShaktiAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    440         /*Default, do nothing*/
     459        /*Update active elements based on ice levelset and ocean levelset*/
     460        GetMaskOfIceVerticesLSMx(femmodel,true);
     461        SetActiveNodesLSMx(femmodel,true);
     462
     463        IssmDouble rho_ice   = femmodel->parameters->FindParam(MaterialsRhoIceEnum);
     464        IssmDouble rho_water = femmodel->parameters->FindParam(MaterialsRhoFreshwaterEnum);
     465        IssmDouble g         = femmodel->parameters->FindParam(ConstantsGEnum);
     466
     467        /*Constrain all nodes that are grounded and unconstrain the ones that float*/
     468        for(Object* & object : femmodel->elements->objects){
     469
     470                /*Get current element and return if not on base*/
     471                Element *element  = xDynamicCast<Element*>(object);
     472                if(!element->IsOnBase()) continue;
     473
     474                int         numnodes  = element->GetNumberOfNodes();
     475                IssmDouble *mask      = xNew<IssmDouble>(numnodes);
     476                IssmDouble *bed       = xNew<IssmDouble>(numnodes);
     477                IssmDouble *thickness = xNew<IssmDouble>(numnodes);
     478                IssmDouble *ls_active = xNew<IssmDouble>(numnodes);
     479
     480                element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum);
     481                element->GetInputListOnNodes(&bed[0],BaseEnum);
     482                element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
     483                element->GetInputListOnNodes(&ls_active[0],HydrologyMaskNodeActivationEnum);
     484
     485                //for(int in=0;in<numnodes;in++){ //
     486                for(int in=0;in<3;in++){ //
     487                        Node* node=element->GetNode(in);
     488                        if(mask[in]>0. && ls_active[in]==1.){
     489                                node->Activate(); //Not sure if we need this!
     490                        }
     491                        else{
     492                                IssmDouble phi =  rho_ice*g*thickness[in] + rho_water*g*bed[in]; //FIXME this is correct!
     493                                node->Deactivate();// Not sure if we need this
     494                                node->ApplyConstraint(0,phi);
     495                        }
     496                }
     497                xDelete<IssmDouble>(mask);
     498                xDelete<IssmDouble>(bed);
     499                xDelete<IssmDouble>(thickness);
     500                xDelete<IssmDouble>(ls_active);
     501        }
     502
    441503        return;
    442504}/*}}}*/
    … …  
    475537
    476538        /*Skip if water or ice shelf element*/
    477         if(element->IsAllFloating()) return;
     539        if(element->IsAllFloating() || !element->IsIceInElement()) return;
     540        if(!element->IsOnBase()) return;
    478541
    479542        /*Intermediaries */
    480         IssmDouble newgap = 0.;
     543        IssmDouble  newgap = 0.;
    481544        IssmDouble  Jdet,meltrate,G,dh[2],B,A,n,dt;
    482545        IssmDouble  gap,bed,thickness,head;
    … …  
    484547        IssmDouble  alpha2,frictionheat;
    485548        IssmDouble* xyz_list = NULL;
    486         IssmDouble  dpressure_water[2],dbed[2],PMPheat,dissipation;
     549        IssmDouble  dpressure_water[2],dbed[2],PMPheat,dissipation;
    487550        IssmDouble q = 0.;
    488         IssmDouble channelization = 0.;
     551        IssmDouble channelization = 0.;
    489552
    490553        /*Retrieve all inputs and parameters*/
    … …  
    551614                if(pressure_water>pressure_ice) pressure_water = pressure_ice;
    552615
    553       /* Compute change in sensible heat due to changes in pressure melting point*/
    554            dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
     616                /* Compute change in sensible heat due to changes in pressure melting point*/
     617                dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
    555618                dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
    556619                dissipation=rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]);
    … …  
    558621                meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
    559622
    560                 element->AddInput(DummyEnum,&meltrate,P0Enum);
    561                 element->AddInput(EsaEmotionEnum,&frictionheat,P0Enum);
    562                 element->AddInput(EsaNmotionEnum,&dissipation,P0Enum);
    563                 element->AddInput(EsaUmotionEnum,&PMPheat,P0Enum);
    564 
     623                element->AddBasalInput(DummyEnum,&meltrate,P0Enum);
     624                element->AddBasalInput(EsaEmotionEnum,&frictionheat,P0Enum);
     625                element->AddBasalInput(EsaNmotionEnum,&dissipation,P0Enum);
     626                element->AddBasalInput(EsaUmotionEnum,&PMPheat,P0Enum);
    565627
    566628                newgap += gauss->weight*Jdet*(gap+dt*(
    567                                         meltrate/rho_ice
    568                                         -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
    569                                         +beta*sqrt(vx*vx+vy*vy)
    570                                         ));
     629                                                meltrate/rho_ice
     630                                                -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
     631                                                +beta*sqrt(vx*vx+vy*vy)
     632                                                ));
    571633
    572634
    … …  
    574636
    575637                /* Compute basal water flux */
    576       q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1]));
     638                q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1]));
    577639
    578640                /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */
    … …  
    590652
    591653        /*Add new gap as an input*/
    592         element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
     654        element->AddBasalInput(HydrologyGapHeightEnum,&newgap,P0Enum);
    593655
    594656        /*Divide by connectivity, add basal flux as an input*/
    595657        q = q/totalweights;
    596         element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum);
     658        element->AddBasalInput(HydrologyBasalFluxEnum,&q,P0Enum);
    597659
    598660        /* Divide by connectivity, add degree of channelization as an input */
    599661        channelization = channelization/totalweights;
    600         element->AddInput(DegreeOfChannelizationEnum,&channelization,P0Enum);
     662        element->AddBasalInput(DegreeOfChannelizationEnum,&channelization,P0Enum);
    601663
    602664        /*Clean up and return*/
    … …  
    616678
    617679        /*Skip if water or ice shelf element*/
    618         if(element->IsAllFloating()) return;
     680        if(element->IsAllFloating() || !element->IsIceInElement()) return;
     681        if(!element->IsOnBase()) return;
    619682
    620683        /*Intermediaries*/
    … …  
    633696        Input* base_input      = element->GetInput(BaseEnum);          _assert_(base_input);
    634697
    635 
    636698   Gauss* gauss=element->NewGauss();
    637699   for (int i=0;i<numnodes;i++){
    … …  
    646708
    647709        /*Add new gap as an input*/
    648         element->AddInput(EffectivePressureEnum,N,element->GetElementType());
     710        element->AddBasalInput(EffectivePressureEnum,N,element->GetElementType());
    649711
    650712        /*Clean up and return*/

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

    r26329 r27649  
    352352
    353353        /*Initialize Load Vector and return if necessary*/
    354         Tria*  tria=(Tria*)element;
     354        Tria*  tria=NULL;
     355        if(element->ObjectEnum()==TriaEnum){
     356                tria = (Tria*)this->element;
     357        }
     358        else if(element->ObjectEnum()==PentaEnum){
     359                tria = (Tria*)this->element->SpawnBasalElement();
     360        }
    355361        _assert_(tria->FiniteElement()==P1Enum);
    356362        if(!tria->IsIceInElement() || tria->IsAllFloating()) return NULL;
    … …  
    380386        /*Clean up and return*/
    381387        delete gauss;
     388        if(tria->IsSpawnedElement()){tria->DeleteMaterials(); delete tria;};
    382389        return pe;
    383390}

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

    r27298 r27649  
    147147                                analysis_enum==HydrologyDCInefficientAnalysisEnum ||
    148148                                analysis_enum==HydrologyDCEfficientAnalysisEnum ||
     149                                analysis_enum==HydrologyShaktiAnalysisEnum ||
     150                                analysis_enum==HydrologyGlaDSAnalysisEnum ||
    149151                                analysis_enum==GLheightadvectionAnalysisEnum ||
    150152                                analysis_enum==LevelsetAnalysisEnum

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

    r27462 r27649  
    6565                /*solid earth considerations:*/
    6666                SolidEarthWaterUpdates(femmodel);
    67 
    6867                delete ug;
    69 
    7068        }
    7169

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

    r27491 r27649  
    117117        if(ishydrology){
    118118                /*We may not be running with ismovingfront so we can't assume LevelsetAnalysis is active*/
    119                 femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
     119                int hydrology_model;
     120                femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
     121                if(hydrology_model==HydrologyshaktiEnum){
     122                        femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum);
     123                }
     124                else if(hydrology_model==HydrologyGlaDSEnum){
     125                        femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
     126                }
     127                else{
     128                        _error_("hydrology model not supported yet");
     129                }
    120130        }else if(isdebris){
    121131                femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);
Note: See TracChangeset for help on using the changeset viewer.