source: issm/oecreview/Archive/22819-23185/ISSM-23034-23035.diff@ 23186

Last change on this file since 23186 was 23186, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22819-23185
File size: 20.7 KB
RevLine 
[23186]1Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 23034)
4+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 23035)
5@@ -65,7 +65,7 @@
6 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
7 int FiniteElement(void);
8 IssmDouble FloatingArea(bool scaled);
9- void FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
10+ void FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure);
11 IssmDouble GetArea3D(void){_error_("not implemented yet!");};
12 IssmDouble GetAreaSpherical(void){_error_("not implemented yet!");};
13 void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
14Index: ../trunk-jpl/src/c/classes/Elements/Element.h
15===================================================================
16--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 23034)
17+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 23035)
18@@ -212,7 +212,7 @@
19 virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
20 virtual int FiniteElement(void)=0;
21 virtual IssmDouble FloatingArea(bool scaled)=0;
22- virtual void FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating)=0;
23+ virtual void FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure)=0;
24 virtual Element* GetBasalElement(void)=0;
25 virtual int GetElementType(void)=0;
26 virtual IssmDouble GetHorizontalSurfaceArea(void){_error_("not implemented");};
27Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
28===================================================================
29--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23034)
30+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23035)
31@@ -1041,7 +1041,7 @@
32 return floatingarea;
33 }
34 /*}}}*/
35-void Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
36+void Tria::FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure){/*{{{*/
37
38 if(!IsOnBase()) return;
39
40@@ -1049,82 +1049,70 @@
41 inputs->GetInputValue(&approximation,ApproximationEnum);
42
43 if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
44- for(int i=0;i<NUMVERTICES;i++){
45- vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
46- vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
47- }
48+ _error_(" contact contiditon only works for FS elements");
49 }
50- else{
51- /*Intermediaries*/
52- IssmDouble* xyz_list = NULL;
53- IssmDouble* xyz_list_base = NULL;
54- IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base;
55- IssmDouble bed_normal[2];
56- IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
57- IssmDouble surface=0,value=0;
58- bool grounded;
59+ /*Intermediaries*/
60+ IssmDouble* xyz_list = NULL;
61+ IssmDouble bed_normal[2],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
62+ IssmDouble water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
63+ IssmDouble sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmaxy[NUMVERTICES],sigma_nn[NUMVERTICES];
64+ IssmDouble viscosity,epsilon[NUMVERTICES];
65+ GetInputListOnVertices(&base[0],BaseEnum);
66+ GetInputListOnVertices(&bed[0],BedEnum);
67+ GetInputListOnVertices(&surface[0],SurfaceEnum);
68+ GetInputListOnVertices(&pressure[0],PressureEnum);
69+ GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
70+ IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
71+ IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
72+ IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum);
73
74- /* Get node coordinates and dof list: */
75- GetVerticesCoordinates(&xyz_list);
76- GetVerticesCoordinatesBase(&xyz_list_base);
77+ /* Get node coordinates and dof list: */
78+ GetVerticesCoordinates(&xyz_list);
79+ /*Retrieve all inputs we will be needing: */
80+ Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
81+ Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
82
83- /*Retrieve all inputs we will be needing: */
84- Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
85- Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
86- Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
87- Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
88- Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
89+ /*1. Recover stresses at the base*/
90+ GaussTria* gauss=new GaussTria();
91+ for (int iv=0;iv<NUMVERTICES;iv++){
92+ gauss->GaussVertex(iv);
93
94- /*Create gauss point in the middle of the basal edge*/
95- Gauss* gauss=NewGaussBase(1);
96- gauss->GaussPoint(0);
97+ /*Compute strain rate viscosity and pressure: */
98+ this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
99+ this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
100+ /*FIXME: this is for Hongju only*/
101+ // pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
102+ // if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv];
103
104- if(!IsFloating()){
105- /*Check for basal force only if grounded and touching GL*/
106- // if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
107- this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
108- this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
109- pressure_input->GetInputValue(&pressure, gauss);
110- base_input->GetInputValue(&base, gauss);
111+ /*Compute Stress*/
112+ sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv];
113+ sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv];
114+ sigmaxy[iv]=2*viscosity*epsilon[2];
115+ }
116
117- /*Compute Stress*/
118- IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
119- IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
120- IssmDouble sigma_xy=2.*viscosity*epsilon[2];
121-
122- /*Get normal vector to the bed */
123- NormalBase(&bed_normal[0],xyz_list_base);
124-
125- /*basalforce*/
126- sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
127-
128- /*Compute water pressure*/
129- IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
130- IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
131- IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum);
132- water_pressure=gravity*rho_water*base;
133-
134- /*Compare basal stress to water pressure and determine whether it should ground*/
135- if (sigma_nn<water_pressure) grounded=true;
136- else grounded=false;
137+ /*2. compute contact condition*/
138+ for(int i=0;i<NUMVERTICES;i++){
139+ /*If was grounded*/
140+ if (phi[i]>=0.){
141+ NormalBase(&bed_normal[0],xyz_list);
142+ sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]);
143+ water_pressure[i]=-gravity*rho_water*base[i];
144+ vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL);
145+ vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL);
146 }
147- else{
148- /*Check for basal elevation if floating*/
149- base_input->GetInputValue(&base, gauss);
150- bed_input->GetInputValue(&bed, gauss);
151- if(base<bed) grounded=true;
152- else grounded=false;
153+ /*If was floating*/
154+ else{
155+ /*Tricky part:
156+ * 1. if base is now touching, we put 1 for sigma_nn and leave water pressure at 0 so that the rest of the module will reground this vertex
157+ * 2. if base is still above bed, water pressure is set as 1, sigma_nn is left as 0, so the GL module will keep it afloat*/
158+ if(base[i]<bed[i]) vertex_sigmann->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
159+ vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
160 }
161- for(int i=0;i<NUMVERTICES;i++){
162- if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
163- else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
164- }
165+ }
166
167- /*clean up*/
168- delete gauss;
169- xDelete<IssmDouble>(xyz_list);
170- xDelete<IssmDouble>(xyz_list_base);
171- }
172+ /*clean up*/
173+ delete gauss;
174+ xDelete<IssmDouble>(xyz_list);
175 }
176 /*}}}*/
177 IssmDouble Tria::GetArea(void){/*{{{*/
178Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
179===================================================================
180--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 23034)
181+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 23035)
182@@ -74,7 +74,7 @@
183 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
184 int FiniteElement(void);
185 IssmDouble FloatingArea(bool scaled);
186- void FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
187+ void FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure);
188 Element* GetBasalElement(void){_error_("not implemented yet");};
189 void GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues);
190 void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
191Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
192===================================================================
193--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23034)
194+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23035)
195@@ -740,90 +740,86 @@
196 return floatingarea;
197 }
198 /*}}}*/
199-void Penta::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
200+void Penta::FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure){/*{{{*/
201
202 if(!IsOnBase()) return;
203
204 int approximation;
205 inputs->GetInputValue(&approximation,ApproximationEnum);
206- if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
207- for(int i=0;i<NUMVERTICES;i++){
208- vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
209- vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
210- }
211+ if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum || approximation==HOFSApproximationEnum){
212+ _error_("Cannot compute contact condition for non FS elements");
213 }
214- else {
215- /*Intermediaries*/
216- IssmDouble* xyz_list = NULL;
217- IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base;
218- IssmDouble bed_normal[3];
219- IssmDouble epsilon[6]; /* epsilon=[exx eyy ezz exy exz eyz];*/
220- IssmDouble surface=0,value=0;
221- bool grounded;
222
223- /* Get node coordinates and dof list: */
224- GetVerticesCoordinates(&xyz_list);
225+ /*Intermediaries*/
226+ IssmDouble* xyz_list = NULL;
227+ IssmDouble bed_normal[3],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
228+ IssmDouble water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
229+ IssmDouble sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmazz[NUMVERTICES],sigmaxy[NUMVERTICES];
230+ IssmDouble sigmayz[NUMVERTICES],sigmaxz[NUMVERTICES],sigma_nn[NUMVERTICES];
231+ IssmDouble viscosity,epsilon[NUMVERTICES];
232+ GetInputListOnVertices(&base[0],BaseEnum);
233+ GetInputListOnVertices(&bed[0],BedEnum);
234+ GetInputListOnVertices(&surface[0],SurfaceEnum);
235+ GetInputListOnVertices(&pressure[0],PressureEnum);
236+ GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
237+ IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
238+ IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
239+ IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum);
240
241- /*Retrieve all inputs we will be needing: */
242- Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
243- Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
244- Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
245- Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
246- Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
247- Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input);
248+ /* Get node coordinates and dof list: */
249+ GetVerticesCoordinates(&xyz_list);
250
251- /*Create gauss point in the middle of the basal edge*/
252- Gauss* gauss=NewGaussBase(1);
253- gauss->GaussPoint(0);
254+ /*Retrieve all inputs we will be needing: */
255+ Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
256+ Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
257+ Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input);
258
259- if(!IsFloating()){
260- /*Check for basal force only if grounded and touching GL*/
261- this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
262- this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
263- pressure_input->GetInputValue(&pressure, gauss);
264- base_input->GetInputValue(&base, gauss); _assert_(base<0.);
265+ /*1. Recover stresses at the base*/
266+ GaussPenta* gauss=new GaussPenta();
267+ for (int iv=0;iv<NUMVERTICES;iv++){
268+ gauss->GaussVertex(iv);
269
270- /*Compute Stress*/
271- IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
272- IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
273- IssmDouble sigma_zz=2.*viscosity*epsilon[2]-pressure;
274- IssmDouble sigma_xy=2.*viscosity*epsilon[3];
275- IssmDouble sigma_xz=2.*viscosity*epsilon[4];
276- IssmDouble sigma_yz=2.*viscosity*epsilon[5];
277+ /*Compute strain rate viscosity and pressure: */
278+ this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
279+ this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
280+ /*FIXME: this is for Hongju only*/
281+ //pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
282+ //if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv];
283
284- /*Get normal vector to the bed */
285- NormalBase(&bed_normal[0],xyz_list);
286+ /*Compute Stress*/
287+ sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv]; // sigma = nu eps - pressure
288+ sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv];
289+ sigmazz[iv]=2*viscosity*epsilon[2]-pressure[iv];
290+ sigmaxy[iv]=2*viscosity*epsilon[3];
291+ sigmaxz[iv]=2*viscosity*epsilon[4];
292+ sigmayz[iv]=2*viscosity*epsilon[5];
293+ }
294
295- /*basalforce*/
296- sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + sigma_zz*bed_normal[2]*bed_normal[2]
297- + 2.*sigma_xy*bed_normal[0]*bed_normal[1] + 2.*sigma_xz*bed_normal[0]*bed_normal[2] + 2.*sigma_yz*bed_normal[1]*bed_normal[2];
298+ /*2. compute contact condition*/
299+ for(int i=0;i<NUMVERTICES;i++){
300
301- /*Compute water pressure*/
302- IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
303- IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
304- IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum);
305- water_pressure=gravity*rho_water*base;
306+ /*If was grounded*/
307+ if (phi[i]>=0.){
308+ NormalBase(&bed_normal[0],xyz_list);
309+ sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1] + sigmazz[i]*bed_normal[2]*bed_normal[2]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]+2*sigmaxz[i]*bed_normal[0]*bed_normal[2]+2*sigmayz[i]*bed_normal[1]*bed_normal[2]);
310+ water_pressure[i]=-gravity*rho_water*base[i];
311+ vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL);
312+ vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL);
313+ }
314
315- /*Compare basal stress to water pressure and determine whether it should ground*/
316- if (sigma_nn<water_pressure) grounded=true;
317- else grounded=false;
318- }
319+ /*If was floating*/
320 else{
321- /*Check for basal elevation if floating*/
322- base_input->GetInputValue(&base, gauss);
323- bed_input->GetInputValue(&bed, gauss);
324- if(base<bed) grounded=true;
325- else grounded=false;
326+ /*Tricky part:
327+ * 1. if base is now touching, we put 1 for sigma_nn and leave water pressure at 0 so that the rest of the module will reground this vertex
328+ * 2. if base is still above bed, water pressure is set as 1, sigma_nn is left as 0, so the GL module will keep it afloat*/
329+ if(base[i]<bed[i]) vertex_sigmann->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
330+ else vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
331 }
332- for(int i=0;i<NUMVERTICES;i++){
333- if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
334- else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
335- }
336+ }
337
338- /*clean up*/
339- delete gauss;
340- xDelete<IssmDouble>(xyz_list);
341- }
342+ /*clean up*/
343+ delete gauss;
344+ xDelete<IssmDouble>(xyz_list);
345 }
346 /*}}}*/
347 void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
348Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
349===================================================================
350--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23034)
351+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23035)
352@@ -2005,8 +2005,8 @@
353 else shelf=false;
354 }
355 else if(migration_style==ContactEnum){
356- if(this->inputs->Max(MaskGroundediceLevelsetEnum) > 0.) shelf=false;
357- else shelf=true;
358+ if(this->inputs->Min(MaskGroundediceLevelsetEnum) < 0.) shelf=true;
359+ else shelf=false;
360 }
361 else if(migration_style==NoneEnum || migration_style==AggressiveMigrationEnum || migration_style==SoftMigrationEnum || migration_style==GroundingOnlyEnum){ //Floating if all nodes are floating
362 if(this->inputs->Min(MaskGroundediceLevelsetEnum) > 0.) shelf=false;
363@@ -2319,7 +2319,7 @@
364 /*go through vertices, and update inputs, considering them to be TriaVertex type: */
365 for(i=0;i<numvertices;i++){
366 /* Contact FS*/
367- if(migration_style == ContactEnum && phi_ungrounding[vertices[i]->Pid()]<10){
368+ if(migration_style == ContactEnum){
369 phi[i]=phi_ungrounding[vertices[i]->Pid()];
370 if(phi[i]>=0.) b[i]=r[i];
371 }
372Index: ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp
373===================================================================
374--- ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23034)
375+++ ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23035)
376@@ -58,42 +58,39 @@
377
378 IssmDouble* ContactFSLevelset(Elements* elements,Vertices* vertices){ /*{{{*/
379
380- Vector<IssmDouble>* vertexgrounded = NULL;
381- Vector<IssmDouble>* vertexfloating = NULL;
382- IssmDouble* serial_vertexgrounded = NULL;
383- IssmDouble* serial_vertexfloating = NULL;
384+ Vector<IssmDouble>* vertex_sigmann = NULL;
385+ Vector<IssmDouble>* vertex_waterpressure = NULL;
386+ IssmDouble* serial_vertex_sigmann = NULL;
387+ IssmDouble* serial_vertex_waterpressure = NULL;
388 IssmDouble* phi = NULL;
389
390 /*Initialize vector with number of vertices*/
391 int numberofvertices = vertices->NumberOfVertices();
392- vertexgrounded = new Vector<IssmDouble>(numberofvertices);
393- vertexfloating = new Vector<IssmDouble>(numberofvertices);
394+ vertex_sigmann = new Vector<IssmDouble>(numberofvertices);
395+ vertex_waterpressure = new Vector<IssmDouble>(numberofvertices);
396 phi = xNew<IssmDouble>(numberofvertices);
397
398 /*Fill vector vertices_potentially_floating: */
399 for(int i=0;i<elements->Size();i++){
400 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
401- element->FSContactMigration(vertexgrounded,vertexfloating);
402+ element->FSContactMigration(vertex_sigmann,vertex_waterpressure);
403 }
404+ /*Assemble vector and serialize */
405+ vertex_sigmann->Assemble();
406+ vertex_waterpressure->Assemble();
407+ serial_vertex_sigmann=vertex_sigmann->ToMPISerial();
408+ serial_vertex_waterpressure=vertex_waterpressure->ToMPISerial();
409
410- /*Assemble vector and serialize */
411- vertexgrounded->Assemble();
412- vertexfloating->Assemble();
413- serial_vertexgrounded=vertexgrounded->ToMPISerial();
414- serial_vertexfloating=vertexfloating->ToMPISerial();
415 for(int i=0;i<numberofvertices;i++){
416- if (serial_vertexgrounded[i]==1. && serial_vertexfloating[i]==1.) phi[i]=0.;
417- else if (serial_vertexgrounded[i]==1) phi[i]=1;
418- else if (serial_vertexfloating[i]==1) phi[i]=-1;
419- else if (serial_vertexgrounded[i]>10) phi[i]=9999;
420- else phi[i]=-9999;
421+ if (serial_vertex_waterpressure[i] > serial_vertex_sigmann[i]) phi[i]=-1;
422+ else phi[i]=1;
423 }
424
425 /*free ressouces and return: */
426- delete vertexgrounded;
427- delete vertexfloating;
428- xDelete<IssmDouble>(serial_vertexgrounded);
429- xDelete<IssmDouble>(serial_vertexfloating);
430+ delete vertex_sigmann;
431+ delete vertex_waterpressure;
432+ xDelete<IssmDouble>(serial_vertex_sigmann);
433+ xDelete<IssmDouble>(serial_vertex_waterpressure);
434
435 return phi;
436 }
Note: See TracBrowser for help on using the repository browser.