Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16611) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16612) @@ -196,4 +196,7 @@ if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum); } + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); + } if(isFS){ iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16611) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16612) @@ -2455,4 +2455,8 @@ input=this->inputs->GetInput(output_enum); break; + case SurfaceNormalVelocityEnum: + this->ComputeSurfaceNormalVelocity(); + input=this->inputs->GetInput(output_enum); + break; default: _error_("input "<SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); + + this->inputs->DeleteInput(SurfaceNormalVelocityEnum); + break; } @@ -2642,4 +2649,51 @@ *(surface_normal+1) = normal[1]/normal_norm; *(surface_normal+2) = normal[2]/normal_norm; +} +/*}}}*/ +/*FUNCTION Tria::ComputeSurfaceNormalVelocity{{{*/ +void Tria::ComputeSurfaceNormalVelocity(){ + + IssmDouble sum,tangential_vector[2],normal_vector[2],time,ux,uy; + IssmDouble normal_velocity[NUMVERTICES],xyz_list[NUMVERTICES][3]; + IssmDouble value[NUMVERTICES],verticesonsurface[NUMVERTICES]; + + for(int iv=0;ivGetInput(SurfaceSlopeXEnum); _assert_(slope_input); + // Input* slope_input= inputs->GetInput(SurfaceEnum); _assert_(slope_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + + /*Get list of nodes on surface*/ + GetInputListOnVertices(&verticesonsurface[0],MeshVertexonsurfaceEnum); + sum = verticesonsurface[0]+verticesonsurface[1]+verticesonsurface[2]; + _assert_(sum==0. || sum==1. || sum==2.); + + /*Compute normal velocity for surface nodes from L2 projected slope*/ + for(int iv=0;ivGaussNode(P1Enum,iv); + slope_input->GetInputValue(&value[iv],gauss); + // slope_input->GetInputDerivativeValue(&value[iv],&xyz_list[0][0],gauss); + vx_input->GetInputValue(&ux,gauss); + vy_input->GetInputValue(&uy,gauss); + tangential_vector[0]=sqrt(1./(pow(value[iv],2.)+1.)); + tangential_vector[1]=value[iv]*tangential_vector[0]; + normal_vector[0]=-1.*tangential_vector[1]; + normal_vector[1]=tangential_vector[0]; + normal_velocity[iv]=ux*normal_vector[0]+uy*normal_vector[1]; + } + } + + delete gauss; + this->inputs->AddInput(new TriaInput(SurfaceNormalVelocityEnum,&normal_velocity[0],P1Enum)); + } /*}}}*/ @@ -3224,5 +3278,5 @@ Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input); smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1 - Total_Smb=rho_ice*base*smb; // smb on element in kg s-1 + Total_Smb=rho_ice*base*smb; // smb on element in kg s-1 /*Return: */ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16611) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16612) @@ -70,4 +70,5 @@ void ComputeStrainRate(Vector* eps); void ComputeStressTensor(); + void ComputeSurfaceNormalVelocity(); void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); @@ -241,5 +242,4 @@ Seg* SpawnSeg(int index1,int index2); void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]); - #ifdef _HAVE_STRESSBALANCE_ ElementMatrix* CreateKMatrixStressbalanceSSA(void); Index: /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16611) +++ /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16612) @@ -26,5 +26,5 @@ case StressbalanceSolutionEnum: - numanalyses=4; + numanalyses=5; analyses=xNew(numanalyses); analyses[0]=StressbalanceAnalysisEnum; @@ -32,4 +32,5 @@ analyses[2]=StressbalanceSIAAnalysisEnum; analyses[3]=L2ProjectionBaseAnalysisEnum; + analyses[4]=ExtrudeFromBaseAnalysisEnum; break; Index: /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 16611) +++ /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 16612) @@ -48,5 +48,5 @@ /*Compute slopes: */ - if(isSIA) surfaceslope_core(femmodel); + if(isSIA || (isFS && meshtype==Mesh2DverticalEnum)) surfaceslope_core(femmodel); if(isFS){ bedslope_core(femmodel); Index: /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16611) +++ /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16612) @@ -32,4 +32,8 @@ solutionsequence_linear(femmodel); } + if(meshtype==Mesh2DverticalEnum){ + femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum); + extrudefrombase_core(femmodel); + } if(save_results){ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16611) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16612) @@ -321,4 +321,5 @@ FreeSurfaceBaseAnalysisEnum, FreeSurfaceTopAnalysisEnum, + SurfaceNormalVelocityEnum, ExtrudeFromBaseAnalysisEnum, ExtrudeFromTopAnalysisEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16611) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16612) @@ -327,4 +327,5 @@ case FreeSurfaceBaseAnalysisEnum : return "FreeSurfaceBaseAnalysis"; case FreeSurfaceTopAnalysisEnum : return "FreeSurfaceTopAnalysis"; + case SurfaceNormalVelocityEnum : return "SurfaceNormalVelocity"; case ExtrudeFromBaseAnalysisEnum : return "ExtrudeFromBaseAnalysis"; case ExtrudeFromTopAnalysisEnum : return "ExtrudeFromTopAnalysis"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16611) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16612) @@ -333,4 +333,5 @@ else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; + else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum; else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; - else if (strcmp(name,"FileParam")==0) return FileParamEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"Input")==0) return InputEnum; + if (strcmp(name,"FileParam")==0) return FileParamEnum; + else if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; - else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressTensor")==0) return StressTensorEnum; + if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; + else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16611) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16612) @@ -319,4 +319,5 @@ def FreeSurfaceBaseAnalysisEnum(): return StringToEnum("FreeSurfaceBaseAnalysis")[0] def FreeSurfaceTopAnalysisEnum(): return StringToEnum("FreeSurfaceTopAnalysis")[0] +def SurfaceNormalVelocityEnum(): return StringToEnum("SurfaceNormalVelocity")[0] def ExtrudeFromBaseAnalysisEnum(): return StringToEnum("ExtrudeFromBaseAnalysis")[0] def ExtrudeFromTopAnalysisEnum(): return StringToEnum("ExtrudeFromTopAnalysis")[0]