Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17207) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17208) @@ -67,19 +67,12 @@ bool save_results; - /* extrapolate */ - Analysis* analysis = new ExtrapolationAnalysis(); - femmodel->parameters->SetParam(VxEnum,ExtrapolationVariableEnum); - analysis->Core(femmodel); - delete analysis; - /*activate formulation: */ femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum); + if(VerboseSolution()) _printf0_("call computational core:\n"); + solutionsequence_linear(femmodel); + /*recover parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - - if(VerboseSolution()) _printf0_("call computational core:\n"); - solutionsequence_linear(femmodel); - if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); @@ -324,2 +317,56 @@ }/*}}}*/ +/* Update of constraints */ +void LevelsetAnalysis::UpdateNoIceConstraints(FemModel* femmodel){/*{{{*/ + + IssmDouble* mask_ice = GetMaskOfIce(femmodel->elements, femmodel->nodes); + + for (int i=0;inodes->Size();i++){ + Node* node=dynamic_cast(femmodel->nodes->GetObjectByOffset(i)); + if(node->InAnalysis(LevelsetAnalysisEnum)){ + if(mask_ice[node->Sid()]==1.){//FIXME: what should be done with actual spcs to ice model? +// node->DofInFSet(0); /*remove spc*/ + } + else{ + IssmDouble defval=0.; + node->ApplyConstraint(1,defval); /*apply spc*/ + } + } + } +}/*}}}*/ +IssmDouble* LevelsetAnalysis::GetMaskOfIce(Elements* elements, Nodes* nodes){/*{{{*/ + + int i; + IssmDouble* mask_ice = NULL; + Vector* vec_mask_ice = NULL; + Element* element = NULL; + + /*Initialize vector with number of vertices*/ + IssmDouble numnodes=nodes->NumberOfNodes(LevelsetAnalysisEnum); + vec_mask_ice=new Vector(numnodes); //nodes at ice front that have ice at next time step + for(i=0;iSize();i++){ + element=dynamic_cast(elements->GetObjectByOffset(i)); + SetMaskOfIceElement(vec_mask_ice, element); + } + + /*Assemble vector and serialize */ + vec_mask_ice->Assemble(); + mask_ice=vec_mask_ice->ToMPISerial(); + delete vec_mask_ice; + return mask_ice; + +}/*}}}*/ +void LevelsetAnalysis::SetMaskOfIceElement(Vector* vec_mask_ice, Element* element){/*{{{*/ + + /* Intermediaries */ + int numnodes = element->GetNumberOfNodes(); + + if(element->IsIceInElement()){ + for(int i = 0;iSetValue(element->nodes[i]->Sid(),1.,INS_VAL); + } + } +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 17207) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 17208) @@ -31,4 +31,8 @@ void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); + /* Updating constraints */ + void UpdateNoIceConstraints(FemModel* femmodel); + IssmDouble* GetMaskOfIce(Elements* elements, Nodes* nodes); + void SetMaskOfIceElement(Vector* vec_mask_ice, Element* element); }; #endif Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17207) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17208) @@ -1368,5 +1368,5 @@ /*If no front, return NULL*/ - if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; + if(!element->IsIcefront()) return NULL; /*Intermediaries*/ @@ -1391,6 +1391,10 @@ IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); element->GetVerticesCoordinates(&xyz_list); - element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + +// element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + element->NormalSection(&normal[0],xyz_list_front); +// element->GetNormalFromLSF(&normal[0]); /*Start looping on Gaussian points*/ @@ -1781,5 +1785,5 @@ /*If no front, return NULL*/ - if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; + if(!element->IsIcefront()) return NULL; /*Intermediaries*/ @@ -2204,5 +2208,5 @@ /*If no front, return NULL*/ - if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; + if(!element->IsIcefront()) return NULL; /*Intermediaries*/ @@ -2970,5 +2974,5 @@ /*If no front, return NULL*/ - if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; + if(!element->IsIcefront()) return NULL; /*Intermediaries*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17207) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17208) @@ -220,4 +220,5 @@ virtual int PressureInterpolation()=0; virtual bool IsZeroLevelset(int levelset_enum)=0; + virtual bool IsIcefront(void)=0; virtual void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0; virtual void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17207) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17208) @@ -114,4 +114,5 @@ int PressureInterpolation(); bool IsZeroLevelset(int levelset_enum); + bool IsIcefront(void){_error_("not implemented yet");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17207) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17208) @@ -159,4 +159,5 @@ void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");}; bool IsZeroLevelset(int levelset_enum){_error_("not implemented");}; + bool IsIcefront(void){_error_("not implemented yet");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");}; void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17207) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17208) @@ -909,12 +909,28 @@ } } - + _assert_(nrfrontnodes==2); + + /* arrange order of frontnodes such that they are oriented counterclockwise */ + if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=2){ + int index=indicesfront[0]; + indicesfront[0]=indicesfront[1]; + indicesfront[1]=index; + } IssmDouble* xyz_front = xNew(3*nrfrontnodes); /* Return nodes */ for(i=0;i(indicesfront); }/*}}}*/ -void Tria::GetNormalFromLSF(IssmDouble *pnormal){/*{{{*/ +void Tria::GetNormalFromLSF(IssmDouble *normal){/*{{{*/ /* Intermediaries */ @@ -929,5 +945,4 @@ IssmDouble* xyz_list = NULL; IssmDouble dlevelset[dim], norm_dlevelset; - IssmDouble normal[dim]={0.}; /*Retrieve all inputs and parameters*/ @@ -936,4 +951,6 @@ counter=0; + for(i=0;iNewGauss(2); for(int ig=gauss->begin();igend();ig++){ @@ -947,7 +964,5 @@ } _assert_(counter>0); - for(i=0;iCore(femmodel); - delete analysis; + + /* extrapolate */ + Analysis* extanalysis = new ExtrapolationAnalysis(); + int vars[2] = {VxEnum, VyEnum}; + for(int iv=0; i<2;i++){ + femmodel->parameters->SetParam(vars[i],ExtrapolationVariableEnum); + extanalysis->Core(femmodel); + } + delete extanalysis; + + LevelsetAnalysis* lsanalysis = new LevelsetAnalysis(); + /* solve level-set equation */ + lsanalysis->Core(femmodel); + /* update spc boundary conditions for ice model */ +// lsanalysis->UpdateNoIceConstraints(femmodel); + delete lsanalysis; }