source: issm/oecreview/Archive/16554-17801/ISSM-16611-16612.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 11.5 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16611)
4+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16612)
5@@ -195,6 +195,9 @@
6 iomodel->FetchDataToInput(elements,VzEnum,0.);
7 if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
8 }
9+ if(iomodel->meshtype==Mesh2DverticalEnum){
10+ iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
11+ }
12 if(isFS){
13 iomodel->FetchDataToInput(elements,MeshVertexonbedEnum);
14 iomodel->FetchDataToInput(elements,PressureEnum,0.);
15Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
16===================================================================
17--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16611)
18+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16612)
19@@ -320,6 +320,7 @@
20 MasstransportSolutionEnum,
21 FreeSurfaceBaseAnalysisEnum,
22 FreeSurfaceTopAnalysisEnum,
23+ SurfaceNormalVelocityEnum,
24 ExtrudeFromBaseAnalysisEnum,
25 ExtrudeFromTopAnalysisEnum,
26 SteadystateSolutionEnum,
27Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
28===================================================================
29--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16611)
30+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16612)
31@@ -326,6 +326,7 @@
32 case MasstransportSolutionEnum : return "MasstransportSolution";
33 case FreeSurfaceBaseAnalysisEnum : return "FreeSurfaceBaseAnalysis";
34 case FreeSurfaceTopAnalysisEnum : return "FreeSurfaceTopAnalysis";
35+ case SurfaceNormalVelocityEnum : return "SurfaceNormalVelocity";
36 case ExtrudeFromBaseAnalysisEnum : return "ExtrudeFromBaseAnalysis";
37 case ExtrudeFromTopAnalysisEnum : return "ExtrudeFromTopAnalysis";
38 case SteadystateSolutionEnum : return "SteadystateSolution";
39Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
40===================================================================
41--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16611)
42+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16612)
43@@ -332,6 +332,7 @@
44 else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
45 else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
46 else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
47+ else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
48 else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
49 else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
50 else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
51@@ -381,11 +382,11 @@
52 else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
53 else if (strcmp(name,"Element")==0) return ElementEnum;
54 else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
55- else if (strcmp(name,"FileParam")==0) return FileParamEnum;
56 else stage=4;
57 }
58 if(stage==4){
59- if (strcmp(name,"Input")==0) return InputEnum;
60+ if (strcmp(name,"FileParam")==0) return FileParamEnum;
61+ else if (strcmp(name,"Input")==0) return InputEnum;
62 else if (strcmp(name,"IntInput")==0) return IntInputEnum;
63 else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum;
64 else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum;
65@@ -504,11 +505,11 @@
66 else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
67 else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
68 else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
69- else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
70 else stage=5;
71 }
72 if(stage==5){
73- if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
74+ if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
75+ else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
76 else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
77 else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
78 else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
79Index: ../trunk-jpl/src/c/cores/surfaceslope_core.cpp
80===================================================================
81--- ../trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16611)
82+++ ../trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16612)
83@@ -31,6 +31,10 @@
84 femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToL2ProjectEnum);
85 solutionsequence_linear(femmodel);
86 }
87+ if(meshtype==Mesh2DverticalEnum){
88+ femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum);
89+ extrudefrombase_core(femmodel);
90+ }
91
92 if(save_results){
93 if(VerboseSolution()) _printf0_("saving results:\n");
94Index: ../trunk-jpl/src/c/cores/stressbalance_core.cpp
95===================================================================
96--- ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 16611)
97+++ ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 16612)
98@@ -47,7 +47,7 @@
99 }
100
101 /*Compute slopes: */
102- if(isSIA) surfaceslope_core(femmodel);
103+ if(isSIA || (isFS && meshtype==Mesh2DverticalEnum)) surfaceslope_core(femmodel);
104 if(isFS){
105 bedslope_core(femmodel);
106 femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
107Index: ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp
108===================================================================
109--- ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16611)
110+++ ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16612)
111@@ -25,12 +25,13 @@
112 switch(solutiontype){
113
114 case StressbalanceSolutionEnum:
115- numanalyses=4;
116+ numanalyses=5;
117 analyses=xNew<int>(numanalyses);
118 analyses[0]=StressbalanceAnalysisEnum;
119 analyses[1]=StressbalanceVerticalAnalysisEnum;
120 analyses[2]=StressbalanceSIAAnalysisEnum;
121 analyses[3]=L2ProjectionBaseAnalysisEnum;
122+ analyses[4]=ExtrudeFromBaseAnalysisEnum;
123 break;
124
125 case SteadystateSolutionEnum:
126Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
127===================================================================
128--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16611)
129+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16612)
130@@ -2454,6 +2454,10 @@
131 this->ComputeStressTensor();
132 input=this->inputs->GetInput(output_enum);
133 break;
134+ case SurfaceNormalVelocityEnum:
135+ this->ComputeSurfaceNormalVelocity();
136+ input=this->inputs->GetInput(output_enum);
137+ break;
138 default:
139 _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
140 }
141@@ -2485,6 +2489,9 @@
142 for(int i=0;i<NUMVERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
143
144 vector->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
145+
146+ this->inputs->DeleteInput(SurfaceNormalVelocityEnum);
147+
148 break;
149 }
150 default:
151@@ -2643,6 +2650,53 @@
152 *(surface_normal+2) = normal[2]/normal_norm;
153 }
154 /*}}}*/
155+/*FUNCTION Tria::ComputeSurfaceNormalVelocity{{{*/
156+void Tria::ComputeSurfaceNormalVelocity(){
157+
158+ IssmDouble sum,tangential_vector[2],normal_vector[2],time,ux,uy;
159+ IssmDouble normal_velocity[NUMVERTICES],xyz_list[NUMVERTICES][3];
160+ IssmDouble value[NUMVERTICES],verticesonsurface[NUMVERTICES];
161+
162+ for(int iv=0;iv<NUMVERTICES;iv++){
163+ normal_velocity[iv]=0.;
164+ value[iv]=0.;
165+ }
166+
167+ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
168+
169+ GaussTria* gauss=new GaussTria();
170+ Input* slope_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slope_input);
171+ // Input* slope_input= inputs->GetInput(SurfaceEnum); _assert_(slope_input);
172+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
173+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
174+
175+
176+ /*Get list of nodes on surface*/
177+ GetInputListOnVertices(&verticesonsurface[0],MeshVertexonsurfaceEnum);
178+ sum = verticesonsurface[0]+verticesonsurface[1]+verticesonsurface[2];
179+ _assert_(sum==0. || sum==1. || sum==2.);
180+
181+ /*Compute normal velocity for surface nodes from L2 projected slope*/
182+ for(int iv=0;iv<NUMVERTICES;iv++){
183+ if(verticesonsurface[iv] == 1){
184+ gauss->GaussNode(P1Enum,iv);
185+ slope_input->GetInputValue(&value[iv],gauss);
186+ // slope_input->GetInputDerivativeValue(&value[iv],&xyz_list[0][0],gauss);
187+ vx_input->GetInputValue(&ux,gauss);
188+ vy_input->GetInputValue(&uy,gauss);
189+ tangential_vector[0]=sqrt(1./(pow(value[iv],2.)+1.));
190+ tangential_vector[1]=value[iv]*tangential_vector[0];
191+ normal_vector[0]=-1.*tangential_vector[1];
192+ normal_vector[1]=tangential_vector[0];
193+ normal_velocity[iv]=ux*normal_vector[0]+uy*normal_vector[1];
194+ }
195+ }
196+
197+ delete gauss;
198+ this->inputs->AddInput(new TriaInput(SurfaceNormalVelocityEnum,&normal_velocity[0],P1Enum));
199+
200+}
201+/*}}}*/
202 /*FUNCTION Tria::TimeAdapt{{{*/
203 IssmDouble Tria::TimeAdapt(void){
204
205@@ -3223,7 +3277,7 @@
206 /*Now get the average SMB over the element*/
207 Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
208 smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1
209- Total_Smb=rho_ice*base*smb; // smb on element in kg s-1
210+ Total_Smb=rho_ice*base*smb; // smb on element in kg s-1
211
212 /*Return: */
213 return Total_Smb;
214Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
215===================================================================
216--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16611)
217+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16612)
218@@ -69,6 +69,7 @@
219 void ComputeBasalStress(Vector<IssmDouble>* sigma_b);
220 void ComputeStrainRate(Vector<IssmDouble>* eps);
221 void ComputeStressTensor();
222+ void ComputeSurfaceNormalVelocity();
223 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
224 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
225 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
226@@ -240,7 +241,6 @@
227 void SetClone(int* minranks);
228 Seg* SpawnSeg(int index1,int index2);
229 void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
230-
231 #ifdef _HAVE_STRESSBALANCE_
232 ElementMatrix* CreateKMatrixStressbalanceSSA(void);
233 ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void);
234Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
235===================================================================
236--- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16611)
237+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16612)
238@@ -318,6 +318,7 @@
239 def MasstransportSolutionEnum(): return StringToEnum("MasstransportSolution")[0]
240 def FreeSurfaceBaseAnalysisEnum(): return StringToEnum("FreeSurfaceBaseAnalysis")[0]
241 def FreeSurfaceTopAnalysisEnum(): return StringToEnum("FreeSurfaceTopAnalysis")[0]
242+def SurfaceNormalVelocityEnum(): return StringToEnum("SurfaceNormalVelocity")[0]
243 def ExtrudeFromBaseAnalysisEnum(): return StringToEnum("ExtrudeFromBaseAnalysis")[0]
244 def ExtrudeFromTopAnalysisEnum(): return StringToEnum("ExtrudeFromTopAnalysis")[0]
245 def SteadystateSolutionEnum(): return StringToEnum("SteadystateSolution")[0]
Note: See TracBrowser for help on using the repository browser.