Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16611) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16612) @@ -195,6 +195,9 @@ iomodel->FetchDataToInput(elements,VzEnum,0.); if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum); } + if(iomodel->meshtype==Mesh2DverticalEnum){ + iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); + } if(isFS){ iomodel->FetchDataToInput(elements,MeshVertexonbedEnum); iomodel->FetchDataToInput(elements,PressureEnum,0.); Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16611) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16612) @@ -320,6 +320,7 @@ MasstransportSolutionEnum, FreeSurfaceBaseAnalysisEnum, FreeSurfaceTopAnalysisEnum, + SurfaceNormalVelocityEnum, ExtrudeFromBaseAnalysisEnum, ExtrudeFromTopAnalysisEnum, SteadystateSolutionEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16611) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16612) @@ -326,6 +326,7 @@ case MasstransportSolutionEnum : return "MasstransportSolution"; case FreeSurfaceBaseAnalysisEnum : return "FreeSurfaceBaseAnalysis"; case FreeSurfaceTopAnalysisEnum : return "FreeSurfaceTopAnalysis"; + case SurfaceNormalVelocityEnum : return "SurfaceNormalVelocity"; case ExtrudeFromBaseAnalysisEnum : return "ExtrudeFromBaseAnalysis"; case ExtrudeFromTopAnalysisEnum : return "ExtrudeFromTopAnalysis"; case SteadystateSolutionEnum : return "SteadystateSolution"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16611) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16612) @@ -332,6 +332,7 @@ else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; 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; else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; 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; else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; 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; else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; Index: ../trunk-jpl/src/c/cores/surfaceslope_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16611) +++ ../trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16612) @@ -31,6 +31,10 @@ femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToL2ProjectEnum); solutionsequence_linear(femmodel); } + if(meshtype==Mesh2DverticalEnum){ + femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum); + extrudefrombase_core(femmodel); + } if(save_results){ if(VerboseSolution()) _printf0_("saving results:\n"); Index: ../trunk-jpl/src/c/cores/stressbalance_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 16611) +++ ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 16612) @@ -47,7 +47,7 @@ } /*Compute slopes: */ - if(isSIA) surfaceslope_core(femmodel); + if(isSIA || (isFS && meshtype==Mesh2DverticalEnum)) surfaceslope_core(femmodel); if(isFS){ bedslope_core(femmodel); femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); Index: ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp =================================================================== --- ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16611) +++ ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16612) @@ -25,12 +25,13 @@ switch(solutiontype){ case StressbalanceSolutionEnum: - numanalyses=4; + numanalyses=5; analyses=xNew(numanalyses); analyses[0]=StressbalanceAnalysisEnum; analyses[1]=StressbalanceVerticalAnalysisEnum; analyses[2]=StressbalanceSIAAnalysisEnum; analyses[3]=L2ProjectionBaseAnalysisEnum; + analyses[4]=ExtrudeFromBaseAnalysisEnum; break; case SteadystateSolutionEnum: Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16611) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16612) @@ -2454,6 +2454,10 @@ this->ComputeStressTensor(); input=this->inputs->GetInput(output_enum); break; + case SurfaceNormalVelocityEnum: + this->ComputeSurfaceNormalVelocity(); + input=this->inputs->GetInput(output_enum); + break; default: _error_("input "<(connectivity[i]); vector->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); + + this->inputs->DeleteInput(SurfaceNormalVelocityEnum); + break; } default: @@ -2643,6 +2650,53 @@ *(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)); + +} +/*}}}*/ /*FUNCTION Tria::TimeAdapt{{{*/ IssmDouble Tria::TimeAdapt(void){ @@ -3223,7 +3277,7 @@ /*Now get the average SMB over the element*/ 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: */ return Total_Smb; Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16611) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16612) @@ -69,6 +69,7 @@ void ComputeBasalStress(Vector* sigma_b); 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); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); @@ -240,7 +241,6 @@ void SetClone(int* minranks); Seg* SpawnSeg(int index1,int index2); void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]); - #ifdef _HAVE_STRESSBALANCE_ ElementMatrix* CreateKMatrixStressbalanceSSA(void); ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void); Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16611) +++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16612) @@ -318,6 +318,7 @@ def MasstransportSolutionEnum(): return StringToEnum("MasstransportSolution")[0] 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] def SteadystateSolutionEnum(): return StringToEnum("SteadystateSolution")[0]