source: issm/oecreview/Archive/14312-15392/ISSM-15222-15223.diff@ 15393

Last change on this file since 15393 was 15393, checked in by Mathieu Morlighem, 12 years ago

NEW: adding Archive/14312-15392 for oecreview
File size: 17.2 KB
RevLine 
[15393]1Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
2===================================================================
3--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15222)
4+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 15223)
5@@ -10,14 +10,14 @@
6 void solutionsequence_hydro_nonlinear(FemModel* femmodel){
7
8 /*solution : */
9- Vector<IssmDouble>* ug=NULL;
10- Vector<IssmDouble>* uf_sed=NULL;
11- Vector<IssmDouble>* uf_epl=NULL;
12+ Vector<IssmDouble>* ug_sed=NULL;
13+ Vector<IssmDouble>* ug_epl=NULL;
14+ Vector<IssmDouble>* uf=NULL;
15 Vector<IssmDouble>* old_uf=NULL;
16- Vector<IssmDouble>* aged_uf_sed=NULL;
17- Vector<IssmDouble>* aged_uf_epl=NULL;
18+ Vector<IssmDouble>* aged_ug_sed=NULL;
19+ Vector<IssmDouble>* aged_ug_epl=NULL;
20 Vector<IssmDouble>* ys=NULL;
21- Vector<IssmDouble>* duf=NULL;
22+ Vector<IssmDouble>* dug=NULL;
23
24 Matrix<IssmDouble>* Kff=NULL;
25 Matrix<IssmDouble>* Kfs=NULL;
26@@ -49,22 +49,21 @@
27
28 /*Iteration on the two layers + transfer*/
29 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
30- GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
31- Reducevectorgtofx(&uf_sed, ug, femmodel->nodes,femmodel->parameters); _assert_(uf_sed);
32+ GetSolutionFromInputsx(&ug_sed, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
33 if(isefficientlayer) {
34 femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
35- GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
36- Reducevectorgtofx(&uf_epl, ug, femmodel->nodes,femmodel->parameters);
37- }
38+ GetSolutionFromInputsx(&ug_epl, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
39+}
40
41 for(;;){
42 sedcount=1;
43 eplcount=1;
44 //save pointer to old velocity
45- delete aged_uf_sed;aged_uf_sed=uf_sed;
46+ delete aged_ug_sed;
47+ aged_ug_sed=ug_sed;
48 if(isefficientlayer){
49- delete aged_uf_epl;
50- aged_uf_epl=uf_epl;
51+ delete aged_ug_epl;
52+ aged_ug_epl=ug_epl;
53 }
54
55 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
56@@ -72,21 +71,21 @@
57 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum);
58 femmodel->UpdateConstraintsx();
59 femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
60- /* femmodel->HydrologyTransferx();
61- */
62+
63 /*Iteration on the sediment layer*/
64 for(;;){
65 femmodel->HydrologyTransferx();
66 femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax);
67 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
68 Reduceloadx(pf,Kfs,ys); delete Kfs;
69- if(sedcount>1)delete uf_sed;
70- Solverx(&uf_sed, Kff, pf,old_uf, df, femmodel->parameters);
71- delete old_uf; old_uf=uf_sed->Duplicate();
72- delete Kff; delete pf; delete ug; delete df;
73+ delete uf;
74+ Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
75+ delete old_uf; old_uf=uf->Duplicate();
76+ if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting aged_ug_sed*/
77+ delete Kff; delete pf; delete df;
78
79- Mergesolutionfromftogx(&ug,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys;
80- InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
81+ Mergesolutionfromftogx(&ug_sed,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
82+ InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_sed);
83 ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
84
85 if (!sedconverged){
86@@ -101,7 +100,7 @@
87 if(sedconverged){
88 femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
89 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sedconverged,ConvergedEnum);
90- InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
91+ InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_sed);
92 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);
93 break;
94 }
95@@ -121,13 +120,14 @@
96 femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL);
97 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
98 Reduceloadx(pf,Kfs,ys); delete Kfs;
99- if(eplcount>1) delete uf_epl;
100- Solverx(&uf_epl, Kff, pf,old_uf, df, femmodel->parameters);
101- delete old_uf; old_uf=uf_epl->Duplicate();
102- delete Kff;delete pf; delete ug; delete df;
103- Mergesolutionfromftogx(&ug,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
104- InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
105-
106+ delete uf;
107+ Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
108+ delete old_uf; old_uf=uf->Duplicate();
109+ if(eplcount>1) delete ug_epl;
110+ delete Kff;delete pf;
111+ delete df;
112+ Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
113+ InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
114 ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
115
116 if (!eplconverged){
117@@ -142,7 +142,7 @@
118 if(eplconverged){
119 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum);
120 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum);
121- InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
122+ InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
123 break;
124 }
125 }
126@@ -150,11 +150,10 @@
127 /*System convergence check*/
128 if(!hydroconverged){
129 //compute norm(du)/norm(u)
130- _assert_(aged_uf_sed); _assert_(uf_sed);
131- duf=uf_sed->Duplicate(); _assert_(duf);
132- aged_uf_sed->Copy(duf);
133- duf->AYPX(uf_sed,-1.0);
134- ndu_sed=duf->Norm(NORM_TWO); nu_sed=aged_uf_sed->Norm(NORM_TWO);
135+ dug=ug_sed->Duplicate(); _assert_(dug);
136+ aged_ug_sed->Copy(dug);
137+ dug->AYPX(ug_sed,-1.0);
138+ ndu_sed=dug->Norm(NORM_TWO); nu_sed=aged_ug_sed->Norm(NORM_TWO);
139 if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
140 if (!xIsNan<IssmDouble>(eps_hyd)){
141 if (!isefficientlayer){
142@@ -168,8 +167,12 @@
143 }
144 }
145 else{
146- duf=aged_uf_epl->Duplicate(); aged_uf_epl->Copy(duf); duf->AYPX(uf_epl,-1.0);
147- ndu_epl=duf->Norm(NORM_TWO); nu_epl=aged_uf_epl->Norm(NORM_TWO);
148+ dug=ug_epl->Duplicate();_assert_(dug);
149+ aged_ug_epl->Copy(dug);_assert_(aged_ug_epl);
150+ dug->AYPX(ug_epl,-1.0);
151+ ndu_epl=dug->Norm(NORM_TWO);
152+ nu_epl=aged_ug_epl->Norm(NORM_TWO);
153+
154 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!");
155 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
156 if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<eps_hyd){
157@@ -192,15 +195,16 @@
158 if(hydroconverged)break;
159 }
160
161- InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
162+ InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_sed);
163+ InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_epl);
164
165 /*Free ressources: */
166- delete ug;
167- delete uf_sed;
168- delete uf_epl;
169+ delete ug_epl;
170+ delete ug_sed;
171+ delete uf;
172 delete old_uf;
173- delete aged_uf_sed;
174- delete aged_uf_epl;
175- delete duf;
176+ delete aged_ug_sed;
177+ delete aged_ug_epl;
178+ delete dug;
179
180 }
181Index: ../trunk-jpl/src/c/classes/DofIndexing.cpp
182===================================================================
183--- ../trunk-jpl/src/c/classes/DofIndexing.cpp (revision 15222)
184+++ ../trunk-jpl/src/c/classes/DofIndexing.cpp (revision 15223)
185@@ -159,20 +159,7 @@
186 else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
187 }
188 /*}}}*/
189-/*FUNCTION DofIndexing::Deactivate{{{*/
190-void DofIndexing::Deactivate(void){
191- this->active = false;
192
193- /*Constrain to 0. at this point*/
194- for(int i=0;i<this->gsize;i++){
195- this->f_set[i] = false;
196- this->s_set[i] = true;
197- this->svalues[i] = 0.;
198- }
199- return;
200-}
201-/*}}}*/
202-
203 /*Some of the Object functionality: */
204 /*FUNCTION DofIndexing::Echo{{{*/
205 void DofIndexing::Echo(void){
206@@ -234,3 +221,29 @@
207 _printf_("\n");
208 }
209 /*}}}*/
210+/*FUNCTION DofIndexing::Deactivate{{{*/
211+void DofIndexing::Deactivate(void){
212+ this->active = false;
213+
214+ /*Constrain to 0. at this point*/
215+ for(int i=0;i<this->gsize;i++){
216+ this->f_set[i] = false;
217+ this->s_set[i] = true;
218+ this->svalues[i] = 0.;
219+ }
220+ return;
221+}
222+/*}}}*/
223+/*FUNCTION DofIndexing::Activate{{{*/
224+void DofIndexing::Activate(void){
225+ this->active = true;
226+
227+ /*Constrain to 0. at this point*/
228+ for(int i=0;i<this->gsize;i++){
229+ this->f_set[i] = true;
230+ this->s_set[i] = false;
231+ this->svalues[i] = 0.;
232+ }
233+ return;
234+}
235+/*}}}*/
236Index: ../trunk-jpl/src/c/classes/FemModel.cpp
237===================================================================
238--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 15222)
239+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 15223)
240@@ -667,8 +667,9 @@
241 /*}}}*/
242 void FemModel::UpdateConstraintsx(void){ /*{{{*/
243
244+ Element *element = NULL;
245 IssmDouble time;
246- int analysis_type;
247+ int analysis_type;
248
249 /*retrieve parameters: */
250 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
251@@ -677,7 +678,13 @@
252 /*start module: */
253 if(VerboseModule()) _printf0_(" Updating constraints for time: " << time << "\n");
254
255- /*First, update dof constraints in nodes, using constraints: */
256+ /*First, Nodes might be activated/deactivated by element*/
257+ for(int i=0;i<elements->Size();i++){
258+ element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
259+ element->UpdateConstraints();
260+ }
261+
262+ /*Second, constraints might be time dependent: */
263 SpcNodesx(nodes,constraints,parameters,analysis_type);
264
265 /*Now, update degrees of freedoms: */
266Index: ../trunk-jpl/src/c/classes/Node.h
267===================================================================
268--- ../trunk-jpl/src/c/classes/Node.h (revision 15222)
269+++ ../trunk-jpl/src/c/classes/Node.h (revision 15223)
270@@ -86,6 +86,7 @@
271 void GetDofList(int* poutdoflist,int approximation_enum,int setenum);
272 void GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum);
273 void FreezeDof(int dof);
274+ void Activate(void);
275 void Deactivate(void);
276 int IsFloating();
277 int IsGrounded();
278Index: ../trunk-jpl/src/c/classes/Elements/Element.h
279===================================================================
280--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 15222)
281+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 15223)
282@@ -134,6 +134,7 @@
283 #ifdef _HAVE_HYDROLOGY_
284 virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
285 virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
286+ virtual void UpdateConstraints(void)=0;
287 #endif
288 };
289 #endif
290Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
291===================================================================
292--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15222)
293+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15223)
294@@ -2819,6 +2819,26 @@
295 this->parameters=NULL;
296 }
297 /*}}}*/
298+/*FUNCTION Tria::UpdateConstraints{{{*/
299+void Tria::UpdateConstraints(void){
300+
301+ int analysis_type;
302+ parameters->FindParam(&analysis_type,AnalysisTypeEnum);
303+
304+ /*Skip if water element*/
305+ if(IsOnWater()) return;
306+
307+ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
308+ switch(analysis_type){
309+ #ifdef _HAVE_HYDROLOGY_
310+ case HydrologyDCEfficientAnalysisEnum:
311+ this->HydrologyUpdateEplDomain();
312+ break;
313+ #endif
314+ }
315+
316+}
317+/*}}}*/
318 /*FUNCTION Tria::UpdatePotentialUngrounding{{{*/
319 int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
320
321@@ -6406,7 +6426,7 @@
322 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
323
324 /*Get the flag to know if the efficient layer is present*/
325- this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
326+ this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
327
328 /*Use the dof list to index into the solution vector: */
329 for(int i=0;i<numdof;i++){
330@@ -6590,6 +6610,30 @@
331 }
332
333 /*}}}*/
334+/*FUNCTION Tria::HydrologyUpdateEplDomain{{{*/
335+void Tria::HydrologyUpdateEplDomain(void){
336+
337+ /*Intermediaries*/
338+ const int numdof = NDOF1 *NUMVERTICES;
339+ bool isefficientlayer;
340+ IssmDouble activeEpl[numdof];
341+
342+ /*Get parameter and Inout list*/
343+ this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
344+ GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
345+
346+ if(isefficientlayer){
347+ for(int i=0;i<numdof;i++){
348+ if(activeEpl[i]==1.0)
349+ nodes[i]->Activate();
350+ else if (activeEpl[i]==0.0)
351+ nodes[i]->Deactivate();
352+ else
353+ _error_("activation value "<< activeEpl[i] <<" not supported");
354+ }
355+ }
356+}
357+/*}}}*/
358 #endif
359
360 #ifdef _HAVE_DAKOTA_
361Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
362===================================================================
363--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 15222)
364+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 15223)
365@@ -249,6 +249,8 @@
366 void GetSolutionFromInputsHydrologyDCInefficient(Vector<IssmDouble>* solution);
367 void GetSolutionFromInputsHydrologyDCEfficient(Vector<IssmDouble>* solution);
368 void CreateHydrologyWaterVelocityInput(void);
369+ void HydrologyUpdateEplDomain(void);
370+ void UpdateConstraints(void);
371 void InputUpdateFromSolutionHydrology(IssmDouble* solution);
372 void InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
373 void InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
374Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
375===================================================================
376--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15222)
377+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15223)
378@@ -3188,6 +3188,13 @@
379 return nflipped;
380 }
381 /*}}}*/
382+/*FUNCTION Penta::UpdateConstraints{{{*/
383+void Penta::UpdateConstraints(void){
384+
385+ /*Do nothing for now*/
386+
387+}
388+/*}}}*/
389 /*FUNCTION Penta::ViscousHeatingCreateInput {{{*/
390 void Penta::ViscousHeatingCreateInput(void){
391
392@@ -9230,7 +9237,6 @@
393 }
394 /*}}}*/
395 #endif
396-
397 #ifdef _HAVE_BALANCED_
398 /*FUNCTION Penta::CreateKMatrixBalancethickness {{{*/
399 ElementMatrix* Penta::CreateKMatrixBalancethickness(void){
400Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
401===================================================================
402--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 15222)
403+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 15223)
404@@ -303,7 +303,7 @@
405 #endif
406
407 #ifdef _HAVE_HYDROLOGY_
408- void CreateHydrologyWaterVelocityInput(void);
409+ void UpdateConstraints(void);
410 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
411 void GetHydrologyTransfer(Vector<IssmDouble>* transfer);
412 #endif
413Index: ../trunk-jpl/src/c/classes/DofIndexing.h
414===================================================================
415--- ../trunk-jpl/src/c/classes/DofIndexing.h (revision 15222)
416+++ ../trunk-jpl/src/c/classes/DofIndexing.h (revision 15223)
417@@ -48,6 +48,7 @@
418 /*}}}*/
419 /*DofIndexing management: {{{*/
420 DofIndexing* Spawn(int* indices, int numindices);
421+ void Activate(void);
422 void Deactivate(void);
423 /*}}}*/
424
425Index: ../trunk-jpl/src/c/classes/Node.cpp
426===================================================================
427--- ../trunk-jpl/src/c/classes/Node.cpp (revision 15222)
428+++ ../trunk-jpl/src/c/classes/Node.cpp (revision 15223)
429@@ -498,6 +498,13 @@
430
431 }
432 /*}}}*/
433+/*FUNCTION Node::Activate{{{*/
434+void Node::Activate(void){
435+
436+ indexing.Activate();
437+
438+}
439+/*}}}*/
440 /*FUNCTION Node::GetApproximation {{{*/
441 int Node::GetApproximation(){
442
Note: See TracBrowser for help on using the repository browser.