Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 27281) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 27282) @@ -233,4 +233,7 @@ }/*}}}*/ ElementMatrix* HydrologyGlaDSAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + /*Skip if water or ice shelf element*/ + if(element->IsAllFloating() || !element->IsIceInElement()) return NULL; /*Intermediaries */ @@ -332,5 +335,5 @@ /*Skip if water or ice shelf element*/ - if(element->IsAllFloating()) return NULL; + if(element->IsAllFloating() || !element->IsIceInElement()) return NULL; /*Intermediaries */ @@ -443,5 +446,45 @@ }/*}}}*/ void HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ - SetActiveNodesLSMx(femmodel); + + /*Update active elements based on ice levelset and ocean levelset*/ + GetMaskOfIceVerticesLSMx(femmodel,true); + SetActiveNodesLSMx(femmodel,true); + + IssmDouble rho_ice = femmodel->parameters->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = femmodel->parameters->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble g = femmodel->parameters->FindParam(ConstantsGEnum); + + /*Constrain all nodes that are grounded and unconstrain the ones that float*/ + for(Object* & object : femmodel->elements->objects){ + Element *element = xDynamicCast(object); + int numnodes = element->GetNumberOfNodes(); + IssmDouble *mask = xNew(numnodes); + IssmDouble *bed = xNew(numnodes); + IssmDouble *thickness = xNew(numnodes); + IssmDouble *ls_active = xNew(numnodes); + + element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum); + element->GetInputListOnNodes(&bed[0],BaseEnum); + element->GetInputListOnNodes(&thickness[0],BaseEnum); + element->GetInputListOnNodes(&ls_active[0],HydrologyMaskNodeActivationEnum); + + for(int in=0;inGetNode(in); + if(mask[in]>0. && ls_active[in]==1.){ + node->Activate(); //Not sure if we need this! + } + else{ + IssmDouble phi = rho_ice*g*thickness[in] + rho_water*g*bed[in]; //FIXME this is correct! + node->Deactivate();// Not sure if we need this + node->ApplyConstraint(0,phi); + } + } + xDelete(mask); + xDelete(bed); + xDelete(thickness); + xDelete(ls_active); + } + + return; }/*}}}*/ @@ -471,7 +514,4 @@ void HydrologyGlaDSAnalysis::UpdateSheetThickness(Element* element){/*{{{*/ - /*Skip if water or ice shelf element*/ - //if(element->IsAllFloating()) return; - /*Intermediaries */ IssmDouble Jdet,vx,vy,ub,h_old,N,h_r,H,b; @@ -484,4 +524,12 @@ /*Initialize new sheet thickness*/ IssmDouble* h_new = xNew(numvertices); + + /*Set to 0 if inactive element*/ + if(element->IsAllFloating() || !element->IsIceInElement()){ + for(int iv=0;ivAddInput(HydrologySheetThicknessEnum,h_new,P1Enum); + xDelete(h_new); + return; + } /*Retrieve all inputs and parameters*/ @@ -545,5 +593,5 @@ /*Force floating ice to have zero sheet thickness*/ - if(element->IsAllFloating()){ + if(element->IsAllFloating() || !element->IsIceInElement()){ for(int iv=0;ivIsAllFloating()) return; /*Intermediary*/ @@ -583,4 +628,12 @@ Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + /*Set to 0 if inactive element*/ + if(element->IsAllFloating() || !element->IsIceInElement()){ + for(int iv=0;ivAddInput(EffectivePressureEnum,N,P1Enum); + xDelete(N); + return; + } + /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(); @@ -607,9 +660,4 @@ if(xIsNan(N[iv])) _error_("NaN found in solution vector"); if(xIsInf(N[iv])) _error_("Inf found in solution vector"); - } - - /*Force floating ice to have zero effective pressure*/ - if(element->IsAllFloating()){ - for(int iv=0;ivGetInput(MaskOceanLevelsetEnum); _assert_(input); Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 27281) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 27282) @@ -603,11 +603,25 @@ /*Initialize Element matrix and return if necessary*/ Tria* tria=(Tria*)element; - if(!tria->IsIceInElement() || tria->IsAllFloating() || this->boundary){ + if(this->boundary){ this->S = 0.; return; } _assert_(tria->FiniteElement()==P1Enum); + + /*Evaluate all fields on center of edge*/ int index1=tria->GetVertexIndex(vertices[0]); int index2=tria->GetVertexIndex(vertices[1]); + GaussTria* gauss=new GaussTria(); + gauss->GaussEdgeCenter(index1,index2); + + /*Set to 0 if inactive*/ + IssmDouble active; + Input* active_input = element->GetInput(HydrologyMaskNodeActivationEnum); _assert_(active_input); + active_input->GetInputValue(&active,gauss); + if(active!=1.){ + this->S = 0.; + delete gauss; + return; + } /*Intermediaries */ @@ -644,8 +658,4 @@ tx = tx/Lt; ty = ty/Lt; - - /*Evaluate fields on center of edge*/ - GaussTria* gauss=new GaussTria(); - gauss->GaussEdgeCenter(index1,index2); /*Get input values at gauss points*/ Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 27281) +++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 27282) @@ -10,10 +10,11 @@ #include "../modules.h" -void SetActiveNodesLSMx(FemModel* femmodel){/*{{{*/ +void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/ /* activate/deactivate nodes for levelset method according to IceMaskNodeActivation */ - /* find parameters */ - int domaintype; - femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); + /*Determine which node activation mask to pull from*/ + int nodeactivationmask = IceMaskNodeActivationEnum; + if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum; + for(Object* & object : femmodel->elements->objects){ @@ -28,6 +29,6 @@ case TaylorHoodEnum:case XTaylorHoodEnum:case LATaylorHoodEnum: case CrouzeixRaviartEnum:case LACrouzeixRaviartEnum:case OneLayerP4zEnum:{ - Input* input=element->GetInput(IceMaskNodeActivationEnum); - if(!input) _error_("Input " << EnumToStringx(IceMaskNodeActivationEnum) << " not found in element"); + Input* input=element->GetInput(nodeactivationmask); + if(!input) _error_("Input " << EnumToStringx(nodeactivationmask) << " not found in element"); /* Start looping on the number of vertices: */ @@ -45,5 +46,5 @@ } default: - element->GetInputListOnNodes(&mask[0],IceMaskNodeActivationEnum); + element->GetInputListOnNodes(&mask[0],nodeactivationmask); break; } @@ -58,5 +59,9 @@ }/*}}}*/ -void GetMaskOfIceVerticesLSMx0(FemModel* femmodel){/*{{{*/ +void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology){/*{{{*/ + + /*Determine which node activation to construct*/ + int nodeactivationmask = IceMaskNodeActivationEnum; + if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum; /*Initialize vector with number of vertices*/ @@ -68,10 +73,23 @@ /*Fill vector with values: */ - for(Object* & object : femmodel->elements->objects){ - Element* element=xDynamicCast(object); - if(element->IsIceInElement()){ - int nbv = element->GetNumberOfVertices(); - for(int iv=0;ivSetValue(element->vertices[iv]->Pid(),1.,INS_VAL); + if(ishydrology){ + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + if(element->IsIceInElement()){ + int nbv = element->GetNumberOfVertices(); + for(int iv=0;ivSetValue(element->vertices[iv]->Pid(),1.,INS_VAL); + } + } + } + } + else{ + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + if(element->IsIceInElement() && element->IsGrounded()){ + int nbv = element->GetNumberOfVertices(); + for(int iv=0;ivSetValue(element->vertices[iv]->Pid(),1.,INS_VAL); + } } } @@ -80,10 +98,21 @@ /*Assemble vector and serialize */ vec_mask_ice->Assemble(); - InputUpdateFromVectorx(femmodel,vec_mask_ice,IceMaskNodeActivationEnum,VertexPIdEnum); + InputUpdateFromVectorx(femmodel,vec_mask_ice,nodeactivationmask,VertexPIdEnum); delete vec_mask_ice; }/*}}}*/ -void GetMaskOfIceVerticesLSMx(FemModel* femmodel){/*{{{*/ +void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/ - femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum); + /*Set configuration to levelset*/ + if(ishydrology){ + /*We may not be running with ismovingfront so we can't assume LevelsetAnalysis is active*/ + femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum); + } + else{ + femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum); + } + + /*Determine which node activation to construct*/ + int nodeactivationmask = IceMaskNodeActivationEnum; + if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum; /*Create vector on gset*/ @@ -97,13 +126,27 @@ Element* element=xDynamicCast(object); - if(element->IsIceInElement()){ - int numnodes = element->GetNumberOfNodes(); - int gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum); - int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum); - IssmDouble* ones = xNew(gsize_local); - for(int n=0;nSetValues(gsize_local,glist_local,ones,INS_VAL); - xDelete(ones); - xDelete(glist_local); + if(ishydrology){ + if(element->IsIceInElement() && element->IsGrounded()){ + int numnodes = element->GetNumberOfNodes(); + int gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum); + int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum); + IssmDouble* ones = xNew(gsize_local); + for(int n=0;nSetValues(gsize_local,glist_local,ones,INS_VAL); + xDelete(ones); + xDelete(glist_local); + } + } + else{ + if(element->IsIceInElement()){ + int numnodes = element->GetNumberOfNodes(); + int gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum); + int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum); + IssmDouble* ones = xNew(gsize_local); + for(int n=0;nSetValues(gsize_local,glist_local,ones,INS_VAL); + xDelete(ones); + xDelete(glist_local); + } } } @@ -120,5 +163,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - element->InputUpdateFromSolutionOneDof(local_ug,IceMaskNodeActivationEnum); + element->InputUpdateFromSolutionOneDof(local_ug,nodeactivationmask); } Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h (revision 27281) +++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h (revision 27282) @@ -8,6 +8,6 @@ #include "../../classes/classes.h" -void SetActiveNodesLSMx(FemModel* femmodel); -void GetMaskOfIceVerticesLSMx0(FemModel* femmodel); -void GetMaskOfIceVerticesLSMx(FemModel* femmodel); +void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology=false); +void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology=false); +void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology=false); #endif /* _SETACTIVENODESLSMX_H*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27281) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27282) @@ -820,4 +820,5 @@ syn keyword cConstant HydrologyWaterVxEnum syn keyword cConstant HydrologyWaterVyEnum +syn keyword cConstant HydrologyMaskNodeActivationEnum syn keyword cConstant IceEnum syn keyword cConstant IceMaskNodeActivationEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27281) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27282) @@ -816,4 +816,5 @@ HydrologyWaterVxEnum, HydrologyWaterVyEnum, + HydrologyMaskNodeActivationEnum, IceEnum, IceMaskNodeActivationEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27281) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27282) @@ -822,4 +822,5 @@ case HydrologyWaterVxEnum : return "HydrologyWaterVx"; case HydrologyWaterVyEnum : return "HydrologyWaterVy"; + case HydrologyMaskNodeActivationEnum : return "HydrologyMaskNodeActivation"; case IceEnum : return "Ice"; case IceMaskNodeActivationEnum : return "IceMaskNodeActivation"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27281) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27282) @@ -813,4 +813,5 @@ syn keyword juliaConstC HydrologyWaterVxEnum syn keyword juliaConstC HydrologyWaterVyEnum +syn keyword juliaConstC HydrologyMaskNodeActivationEnum syn keyword juliaConstC IceEnum syn keyword juliaConstC IceMaskNodeActivationEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27281) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27282) @@ -840,4 +840,5 @@ else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; + else if (strcmp(name,"HydrologyMaskNodeActivation")==0) return HydrologyMaskNodeActivationEnum; else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Misfit")==0) return MisfitEnum; else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum; - else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; + if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum; + else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; else if (strcmp(name,"NewDamage")==0) return NewDamageEnum; else if (strcmp(name,"Node")==0) return NodeEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum; else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum; - else if (strcmp(name,"SmbD")==0) return SmbDEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; + if (strcmp(name,"SmbD")==0) return SmbDEnum; + else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum; else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"Surface")==0) return SurfaceEnum; else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; - else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; - else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; + if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; + else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum; @@ -1366,9 +1367,9 @@ else if (strcmp(name,"ElementHook")==0) return ElementHookEnum; else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; - else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum; + if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; + else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum; else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum; else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum; @@ -1489,9 +1490,9 @@ else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; - else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; + if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; + else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; else if (strcmp(name,"Materials")==0) return MaterialsEnum; @@ -1612,9 +1613,9 @@ else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum; - else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; + if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum; + else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; else if (strcmp(name,"Sset")==0) return SsetEnum;