source: issm/oecreview/Archive/21724-22754/ISSM-21930-21931.diff

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

CHG: added 21724-22754
File size: 32.2 KB
RevLine 
[22755]1Index: ../trunk-jpl/src/c/bamg/GeomEdge.h
2===================================================================
3--- ../trunk-jpl/src/c/bamg/GeomEdge.h (revision 21930)
4+++ ../trunk-jpl/src/c/bamg/GeomEdge.h (revision 21931)
5@@ -26,7 +26,6 @@
6
7 //Methods
8 R2 F(double theta) const ; // parametrization of the curve edge
9- double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
10 int Cracked() const;
11 int TgA() const;
12 int TgB() const;
13Index: ../trunk-jpl/src/c/bamg/GeomEdge.cpp
14===================================================================
15--- ../trunk-jpl/src/c/bamg/GeomEdge.cpp (revision 21930)
16+++ ../trunk-jpl/src/c/bamg/GeomEdge.cpp (revision 21931)
17@@ -62,77 +62,6 @@
18 int GeomEdge::Mark() const {/*{{{*/
19 return type &16;
20 }/*}}}*/
21- double GeomEdge::R1tg(double theta,R2 & t) const{/*{{{*/
22- /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/R1tg)*/
23- // 1/R of radius of cuvature
24-
25- R2 A=v[0]->r,B=v[1]->r;
26- double dca,dcb,dcta,dctb;
27- double ddca,ddcb,ddcta,ddctb;
28- double tt = theta*theta;
29-
30- //check theta
31- _assert_(theta>=0 && theta<=1);
32-
33- if (TgA()){
34- if (TgB()){
35- // Tangent A and B provided:
36- // interpolation d'hermite
37- dcb = 6*theta*(1-theta);
38- ddcb = 6*(1-2*theta);
39- dca = -dcb;
40- ddca = -ddcb;
41- dcta = (3*theta - 4)*theta + 1;
42- ddcta=6*theta-4;
43- dctb = 3*tt - 2*theta;
44- ddctb = 6*theta-2;
45- }
46- else {
47- //Tangent A provided but tangent B not provided
48- // 1-t*t, t-t*t, t*t
49- double t = theta;
50- dcb = 2*t;
51- ddcb = 2;
52- dca = -dcb;
53- ddca = -2;
54- dcta = 1-dcb;
55- ddcta = -ddcb;
56- dctb=0;
57- ddctb=0;
58- }
59- }
60- else{
61- if (TgB()){
62- //Tangent B provided but tangent A not provided
63- double t = 1-theta;
64- dca = -2*t;
65- ddca = 2;
66- dcb = -dca;
67- ddcb = -2;
68- dctb = 1+dca;
69- ddctb= ddca;
70- dcta =0;
71- ddcta =0;
72- }
73- else {
74- //Neither thangent A nor tangent B provided
75- // lagrange P1
76- t=B-A;
77- return 0;
78- }
79- }
80- R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb;
81- R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb;
82- double d2=(d,d);
83- double sd2 = sqrt(d2);
84- t=d;
85- if(d2>1.0e-20){
86- t/=sd2;
87- return Abs(Det(d,dd))/(d2*sd2);
88- }
89- else return 0;
90- }
91- /*}}}*/
92 int GeomEdge::Required() {/*{{{*/
93 return type &64;
94 }/*}}}*/
95Index: ../trunk-jpl/src/c/bamg/Mesh.cpp
96===================================================================
97--- ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 21930)
98+++ ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 21931)
99@@ -1593,7 +1593,10 @@
100 long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]);
101 long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
102 k = edge4->SortAndFind(i0,i1);
103- _assert_(k>=0);
104+ if(k<0){
105+ delete [] colorV;
106+ _error_("This should not happen");
107+ }
108 subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
109 subdomains[i].edge = edges+k;
110 Gh.subdomains[i].edge = Gh.edges + k;
111@@ -2053,7 +2056,15 @@
112 if (xy==0) dd=dxdx;
113 else if (xy==1) dd=dxdy;
114 else if (xy==2) dd=dydy;
115- else _error_("not supported yet");
116+ else{
117+ delete [] detT;
118+ delete [] Mmass;
119+ delete [] workT;
120+ delete [] workV;
121+ delete [] Mmassxx;
122+ delete [] OnBoundary;
123+ _error_("not supported yet");
124+ }
125 // do leat 2 iteration for boundary problem
126 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
127 for (i=0;i<nbt;i++)
128@@ -2102,10 +2113,10 @@
129 delete [] dxdx;
130 delete [] dxdy;
131 delete [] dydy;
132- delete [] workT;
133+ delete [] workT;
134 delete [] workV;
135 delete [] Mmassxx;
136- delete [] OnBoundary;
137+ delete [] OnBoundary;
138
139 }
140 /*}}}*/
141@@ -2393,6 +2404,7 @@
142 it++;} // end while (it<nbt)
143 if (nbt == nbtout || !NbSubDomTot) {
144 delete [] HeapArete;
145+ delete [] HeapTriangle;
146 _error_("The boundary is not close: all triangles are outside");
147 }
148
149@@ -2882,7 +2894,10 @@
150 if(kk<10) _printf_("BUG: the geometrical edge " << i << " is on no edge curve\n");
151 }
152 }
153- if(kk) _error_("See above");
154+ if(kk){
155+ delete [] e;
156+ _error_("See above");
157+ }
158
159 return e;
160 }
161@@ -4397,6 +4412,7 @@
162 for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
163 printf("\n");
164 if(NbVerticesOnGeomVertex >= maxnbv){
165+ delete [] bcurve;
166 _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
167 }
168
169Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
170===================================================================
171--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21930)
172+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21931)
173@@ -984,7 +984,6 @@
174 GiaIvinsAnalysisEnum,
175 EsaSolutionEnum,
176 EsaAnalysisEnum,
177- MeshdeformationAnalysisEnum,
178 LevelsetAnalysisEnum,
179 ExtrapolationAnalysisEnum,
180 /*}}}*/
181Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
182===================================================================
183--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21930)
184+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21931)
185@@ -952,7 +952,6 @@
186 case GiaIvinsAnalysisEnum : return "GiaIvinsAnalysis";
187 case EsaSolutionEnum : return "EsaSolution";
188 case EsaAnalysisEnum : return "EsaAnalysis";
189- case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis";
190 case LevelsetAnalysisEnum : return "LevelsetAnalysis";
191 case ExtrapolationAnalysisEnum : return "ExtrapolationAnalysis";
192 case ApproximationEnum : return "Approximation";
193Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
194===================================================================
195--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21930)
196+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21931)
197@@ -973,7 +973,6 @@
198 else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum;
199 else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum;
200 else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
201- else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
202 else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
203 else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
204 else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
205@@ -997,11 +996,11 @@
206 else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
207 else if (strcmp(name,"MinVel")==0) return MinVelEnum;
208 else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
209+ else if (strcmp(name,"MinVx")==0) return MinVxEnum;
210 else stage=9;
211 }
212 if(stage==9){
213- if (strcmp(name,"MinVx")==0) return MinVxEnum;
214- else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
215+ if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
216 else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
217 else if (strcmp(name,"MinVy")==0) return MinVyEnum;
218 else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
219Index: ../trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp
220===================================================================
221--- ../trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp (revision 21930)
222+++ ../trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp (revision 21931)
223@@ -1690,20 +1690,3 @@
224 xDelete<IssmPDouble>(work);
225
226 }/*}}}*/
227-
228-/*Element Gauss points TO BE REMOVED*/
229-void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus ) {/*{{{*/
230- /*Gauss quadrature points for the quadrilaterial.*/
231-
232- /*get the gauss points using the product of two line rules */
233- GaussLegendreLinear(pxgaus, pxwgt, nigaus);
234- GaussLegendreLinear(pegaus, pewgt, njgaus);
235-}/*}}}*/
236-void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus ) {/*{{{*/
237- /*Gauss quadrature points for the hexahedron.*/
238-
239- /* get the gauss points using the product of three line rules */
240- GaussLegendreLinear(pxgaus, pxwgt, nigaus);
241- GaussLegendreLinear(pegaus, pewgt, njgaus);
242- GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
243-}/*}}}*/
244Index: ../trunk-jpl/src/c/shared/Numerics/GaussPoints.h
245===================================================================
246--- ../trunk-jpl/src/c/shared/Numerics/GaussPoints.h (revision 21930)
247+++ ../trunk-jpl/src/c/shared/Numerics/GaussPoints.h (revision 21931)
248@@ -17,7 +17,4 @@
249 #define MAX_GAUS_ITER 30
250 void GaussRecur(IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta);
251
252-void gaussQuad(IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus);
253-void gaussHexa(IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus);
254-
255 #endif
256Index: ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp
257===================================================================
258--- ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 21930)
259+++ ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 21931)
260@@ -67,18 +67,8 @@
261 }
262
263 }
264-double ddistance(double x1,double y1,double x2,double y2){
265- return sqrt(pow(x2-x1,2)+pow(y2-y1,2));
266-}
267-double ddot(double x1, double y1, double x2, double y2){
268- return x1*x2+y1*y2;
269-}
270-
271-bool isPointLeftOfRay(double x, double y, double raySx, double raySy, double rayEx, double rayEy) {
272- return (y-raySy)*(rayEx-raySx)
273- > (x-raySx)*(rayEy-raySy);
274-}
275-
276+double ddistance(double x1,double y1,double x2,double y2){ return sqrt(pow(x2-x1,2)+pow(y2-y1,2)); }
277+double ddot(double x1, double y1, double x2, double y2){ return x1*x2+y1*y2; }
278 double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0){
279
280 // Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)
281Index: ../trunk-jpl/src/c/Makefile.am
282===================================================================
283--- ../trunk-jpl/src/c/Makefile.am (revision 21930)
284+++ ../trunk-jpl/src/c/Makefile.am (revision 21931)
285@@ -293,7 +293,6 @@
286 ./cores/dummy_core.cpp\
287 ./cores/surfaceslope_core.cpp\
288 ./cores/bedslope_core.cpp\
289- ./cores/meshdeformation_core.cpp\
290 ./cores/damage_core.cpp\
291 ./cores/levelsetfunctionslope_core.cpp\
292 ./cores/movingfront_core.cpp\
293@@ -447,9 +446,6 @@
294 if SMOOTH
295 issm_sources += ./analyses/SmoothAnalysis.cpp
296 endif
297-if MESHDEFORMATION
298-issm_sources += ./analyses/MeshdeformationAnalysis.cpp
299-endif
300 if LEVELSET
301 issm_sources += ./analyses/LevelsetAnalysis.cpp
302 issm_sources += ./modules/Calvingx/Calvingx.cpp
303Index: ../trunk-jpl/src/c/cores/meshdeformation_core.cpp
304===================================================================
305--- ../trunk-jpl/src/c/cores/meshdeformation_core.cpp (revision 21930)
306+++ ../trunk-jpl/src/c/cores/meshdeformation_core.cpp (nonexistent)
307@@ -1,30 +0,0 @@
308-/*!\file: meshdeformation_core.cpp
309- * \brief: core of the slope solution
310- */
311-
312-#include "./cores.h"
313-#include "../toolkits/toolkits.h"
314-#include "../classes/classes.h"
315-#include "../shared/shared.h"
316-#include "../solutionsequences/solutionsequences.h"
317-#include "../modules/modules.h"
318-
319-void meshdeformation_core(FemModel* femmodel){
320-
321- /*parameters: */
322- bool save_results;
323-
324- /*Recover some parameters: */
325- femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
326-
327- if(VerboseSolution()) _printf0_("computing mesh deformation (elasticity)...\n");
328-
329- /*Call on core computations: */
330- femmodel->SetCurrentConfiguration(MeshdeformationAnalysisEnum);
331- solutionsequence_linear(femmodel);
332-
333- if(save_results){
334- if(VerboseSolution()) _printf0_("saving results:\n");
335- }
336-
337-}
338Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
339===================================================================
340--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21930)
341+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21931)
342@@ -1221,49 +1221,6 @@
343 return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
344 }
345 /*}}}*/
346-void Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
347-
348- /*Intermediary*/
349- int num_controls;
350- int* control_type=NULL;
351- Input* input=NULL;
352-
353- /*retrieve some parameters: */
354- this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
355- this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
356-
357- for(int i=0;i<num_controls;i++){
358-
359- if(control_type[i]==MaterialsRheologyBbarEnum){
360- if (!IsOnBase()) goto cleanup_and_return;
361- input=(Input*)this->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
362- }
363- else if(control_type[i]==DamageDbarEnum){
364- if (!IsOnBase()) goto cleanup_and_return;
365- input=(Input*)this->inputs->GetInput(DamageDEnum); _assert_(input);
366- }
367- else{
368- input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
369- }
370- if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
371-
372- ((ControlInput*)input)->UpdateValue(scalar);
373- ((ControlInput*)input)->Constrain();
374- if (save_parameter) ((ControlInput*)input)->SaveValue();
375-
376- if(control_type[i]==MaterialsRheologyBbarEnum){
377- this->InputExtrude(MaterialsRheologyBEnum,-1);
378- }
379- else if(control_type[i]==DamageDbarEnum){
380- this->InputExtrude(DamageDEnum,-1);
381- }
382- }
383-
384- /*Clean up and return*/
385-cleanup_and_return:
386- xDelete<int>(control_type);
387-}
388-/*}}}*/
389 void Penta::InputDepthAverageAtBase(int original_enum,int average_enum){/*{{{*/
390
391 IssmDouble Jdet,value;
392@@ -1367,18 +1324,6 @@
393 }
394 }
395 /*}}}*/
396-void Penta::InputScale(int enum_type,IssmDouble scale_factor){/*{{{*/
397-
398- Input* input=NULL;
399-
400- /*Make a copy of the original input: */
401- input=(Input*)this->inputs->GetInput(enum_type);
402- if(!input)_error_("could not find old input with enum: " << EnumToStringx(enum_type));
403-
404- /*Scale: */
405- input->Scale(scale_factor);
406-}
407-/*}}}*/
408 void Penta::InputUpdateFromIoModel(int index,IoModel* iomodel){ /*{{{*/
409
410 /*Intermediaries*/
411Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
412===================================================================
413--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 21930)
414+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 21931)
415@@ -90,10 +90,8 @@
416 IssmDouble GroundedArea(void);
417 IssmDouble IceVolume(void);
418 IssmDouble IceVolumeAboveFloatation(void);
419- void InputControlUpdate(IssmDouble scalar,bool save_parameter);
420 void InputDepthAverageAtBase(int enum_type,int average_enum_type);
421 void InputExtrude(int enum_type,int start);
422- void InputScale(int enum_type,IssmDouble scale_factor);
423 void InputUpdateFromIoModel(int index, IoModel* iomodel);
424 void InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
425 void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);
426Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
427===================================================================
428--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 21930)
429+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 21931)
430@@ -78,10 +78,8 @@
431 IssmDouble GroundedArea(void){_error_("not implemented yet");};
432 IssmDouble IceVolume(void){_error_("not implemented yet");};
433 IssmDouble IceVolumeAboveFloatation(void){_error_("not implemented yet");};
434- void InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
435 void InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};
436 void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
437- void InputScale(int enum_type,IssmDouble scale_factor){_error_("not implemented yet");};
438 void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
439 void InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
440 void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
441Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
442===================================================================
443--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21930)
444+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21931)
445@@ -92,10 +92,8 @@
446 bool IsOnBase();
447 bool IsOnSurface();
448 bool IsZeroLevelset(int levelset_enum){_error_("not implemented");};
449- void InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
450 void InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};
451 void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
452- void InputScale(int enum_type,IssmDouble scale_factor){_error_("not implemented yet");};
453 void InputUpdateFromIoModel(int index, IoModel* iomodel);
454 void InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum);
455 void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
456Index: ../trunk-jpl/src/c/classes/Elements/Element.h
457===================================================================
458--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 21930)
459+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 21931)
460@@ -221,10 +221,8 @@
461 virtual IssmDouble GroundedArea(void)=0;
462 virtual IssmDouble IceVolume(void)=0;
463 virtual IssmDouble IceVolumeAboveFloatation(void)=0;
464- virtual void InputControlUpdate(IssmDouble scalar,bool save_parameter)=0;
465 virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
466 virtual void InputExtrude(int input_enum,int start)=0;
467- virtual void InputScale(int enum_type,IssmDouble scale_factor)=0;
468 virtual void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0;
469 virtual void InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum)=0;
470 virtual bool IsFaceOnBoundary(void)=0;
471Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
472===================================================================
473--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21930)
474+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21931)
475@@ -1366,7 +1366,7 @@
476 /*Computeportion of the element that has a positive levelset*/
477
478 bool negative=true;
479- int point;
480+ int point=0;
481 const IssmPDouble epsilon= 1.e-15;
482 IssmDouble f1,f2;
483
484@@ -1404,6 +1404,9 @@
485 f1=gl[1]/(gl[1]-gl[2]);
486 f2=gl[1]/(gl[1]-gl[0]);
487 }
488+ else{
489+ _error_("This case should NOT be happening");
490+ }
491 }
492 *point1=point;
493 *fraction1=f1;
494@@ -1660,33 +1663,6 @@
495 return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
496 }
497 /*}}}*/
498-void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
499-
500- /*Intermediary*/
501- int num_controls;
502- int* control_type=NULL;
503- Input* input=NULL;
504-
505- /*retrieve some parameters: */
506- this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
507- this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
508-
509- for(int i=0;i<num_controls;i++){
510- input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
511- if (input->ObjectEnum()!=ControlInputEnum){
512- _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
513- }
514-
515- ((ControlInput*)input)->UpdateValue(scalar);
516- ((ControlInput*)input)->Constrain();
517- if (save_parameter) ((ControlInput*)input)->SaveValue();
518-
519- }
520-
521- /*Clean up and return*/
522- xDelete<int>(control_type);
523-}
524-/*}}}*/
525 void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
526
527 /*New input*/
528@@ -1705,18 +1681,6 @@
529 this->inputs->AddInput((Input*)newinput);
530 }
531 /*}}}*/
532-void Tria::InputScale(int enum_type,IssmDouble scale_factor){/*{{{*/
533-
534- Input* input=NULL;
535-
536- /*Make a copy of the original input: */
537- input=(Input*)this->inputs->GetInput(enum_type);
538- if(!input)_error_("could not find old input with enum: " << EnumToStringx(enum_type));
539-
540- /*Scale: */
541- input->Scale(scale_factor);
542-}
543-/*}}}*/
544 void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index/*{{{*/
545
546 /*Intermediaries*/
547Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
548===================================================================
549--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 21930)
550+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 21931)
551@@ -90,10 +90,8 @@
552 bool HasEdgeOnSurface();
553 IssmDouble IceVolume(void);
554 IssmDouble IceVolumeAboveFloatation(void);
555- void InputControlUpdate(IssmDouble scalar,bool save_parameter);
556 void InputDepthAverageAtBase(int enum_type,int average_enum_type);
557 void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
558- void InputScale(int enum_type,IssmDouble scale_factor);
559 bool IsFaceOnBoundary(void);
560 bool IsIcefront(void);
561 bool IsNodeOnShelfFromFlags(IssmDouble* flags);
562Index: ../trunk-jpl/src/c/classes/Inputs/ControlInput.h
563===================================================================
564--- ../trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 21930)
565+++ ../trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 21931)
566@@ -80,7 +80,6 @@
567 void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
568 void SaveValue(void);
569 void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
570- void ScaleGradient(IssmDouble scale);
571 void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
572 void SetGradient(Input* gradient_in);
573 void SetInput(Input* in_input);
574Index: ../trunk-jpl/src/c/classes/Inputs/DatasetInput.h
575===================================================================
576--- ../trunk-jpl/src/c/classes/Inputs/DatasetInput.h (revision 21930)
577+++ ../trunk-jpl/src/c/classes/Inputs/DatasetInput.h (revision 21931)
578@@ -76,7 +76,6 @@
579 void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
580 void SaveValue(void){_error_("not implemented yet");};
581 void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
582- void ScaleGradient(IssmDouble scale){_error_("not implemented yet");};
583 void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
584 void SetGradient(Input* gradient_in){_error_("not implemented yet");};
585 void UpdateValue(IssmDouble scalar){_error_("not implemented yet");};
586Index: ../trunk-jpl/src/c/classes/Inputs/ControlInput.cpp
587===================================================================
588--- ../trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 21930)
589+++ ../trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 21931)
590@@ -235,10 +235,6 @@
591 if(savedvalues) delete this->savedvalues;
592 this->savedvalues=xDynamicCast<Input*>(this->values->copy());
593 }/*}}}*/
594-void ControlInput::ScaleGradient(IssmDouble scaling_factor){/*{{{*/
595- if(!gradient) _error_("Gradient of ControlInput " << EnumToStringx(enum_type) << " not found");
596- gradient->Scale(scaling_factor);
597-}/*}}}*/
598 void ControlInput::SetGradient(Input* gradient_in){/*{{{*/
599
600 /*Get enum for current gradient*/
601Index: ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp
602===================================================================
603--- ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp (revision 21930)
604+++ ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp (nonexistent)
605@@ -1,56 +0,0 @@
606-#include "./MeshdeformationAnalysis.h"
607-#include "../toolkits/toolkits.h"
608-#include "../classes/classes.h"
609-#include "../shared/shared.h"
610-#include "../modules/modules.h"
611-
612-/*Model processing*/
613-void MeshdeformationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
614- _error_("not implemented yet");
615-}/*}}}*/
616-void MeshdeformationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
617- _error_("not implemented yet");
618-}/*}}}*/
619-void MeshdeformationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
620- _error_("not implemented yet");
621-}/*}}}*/
622-int MeshdeformationAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
623- _error_("not implemented");
624-}/*}}}*/
625-void MeshdeformationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
626- _error_("not implemented yet");
627-}/*}}}*/
628-void MeshdeformationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
629- _error_("not implemented yet");
630-}/*}}}*/
631-
632-/*Finite Element Analysis*/
633-void MeshdeformationAnalysis::Core(FemModel* femmodel){/*{{{*/
634- _error_("not implemented");
635-}/*}}}*/
636-ElementVector* MeshdeformationAnalysis::CreateDVector(Element* element){/*{{{*/
637- /*Default, return NULL*/
638- return NULL;
639-}/*}}}*/
640-ElementMatrix* MeshdeformationAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
641-_error_("Not implemented");
642-}/*}}}*/
643-ElementMatrix* MeshdeformationAnalysis::CreateKMatrix(Element* element){/*{{{*/
644- _error_("not implemented yet");
645-}/*}}}*/
646-ElementVector* MeshdeformationAnalysis::CreatePVector(Element* element){/*{{{*/
647-_error_("not implemented yet");
648-}/*}}}*/
649-void MeshdeformationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
650- _error_("not implemented yet");
651-}/*}}}*/
652-void MeshdeformationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
653- _error_("Not implemented yet");
654-}/*}}}*/
655-void MeshdeformationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
656- _error_("not implemented yet");
657-}/*}}}*/
658-void MeshdeformationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
659- /*Default, do nothing*/
660- return;
661-}/*}}}*/
662Index: ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h
663===================================================================
664--- ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h (revision 21930)
665+++ ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h (nonexistent)
666@@ -1,33 +0,0 @@
667-/*! \file MeshdeformationAnalysis.h
668- * \brief: header file for generic external result object
669- */
670-
671-#ifndef _MeshdeformationAnalysis_
672-#define _MeshdeformationAnalysis_
673-
674-/*Headers*/
675-#include "./Analysis.h"
676-
677-class MeshdeformationAnalysis: public Analysis{
678-
679- public:
680- /*Model processing*/
681- void CreateConstraints(Constraints* constraints,IoModel* iomodel);
682- void CreateLoads(Loads* loads, IoModel* iomodel);
683- void CreateNodes(Nodes* nodes,IoModel* iomodel);
684- int DofsPerNode(int** doflist,int domaintype,int approximation);
685- void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
686- void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
687-
688- /*Finite element Analysis*/
689- void Core(FemModel* femmodel);
690- ElementVector* CreateDVector(Element* element);
691- ElementMatrix* CreateJacobianMatrix(Element* element);
692- ElementMatrix* CreateKMatrix(Element* element);
693- ElementVector* CreatePVector(Element* element);
694- void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
695- void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
696- void InputUpdateFromSolution(IssmDouble* solution,Element* element);
697- void UpdateConstraints(FemModel* femmodel);
698-};
699-#endif
700Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
701===================================================================
702--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 21930)
703+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 21931)
704@@ -551,62 +551,6 @@
705 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
706 }
707 }/*}}}*/
708-void LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
709-
710- /* Intermediaries */
711- int i,k;
712- IssmDouble dmaxp=0.,dmaxm=0,val=0.;
713-
714- /*Initialize vector with number of vertices*/
715- int numvertices=femmodel->vertices->NumberOfVertices();
716- Element* element;
717-
718- Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
719- GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
720-
721- /* set NaN on elements intersected by zero levelset */
722- for(i=0;i<femmodel->elements->Size();i++){
723- element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
724- if(element->IsZeroLevelset(MaskIceLevelsetEnum))
725- for(k=0;k<element->GetNumberOfVertices();k++)
726- vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL);
727- }
728-
729- /* set distance on elements intersected by zero levelset */
730- for(i=0;i<femmodel->elements->Size();i++){
731- element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
732- if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
733- SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
734-
735- /* Get maximum distance to interface along vertices */
736- for(k=0;k<element->GetNumberOfVertices();k++){
737- vec_dist_zerolevelset->GetValue(&val,element->vertices[k]->Sid());
738- if((val>0.) && (val>dmaxp))
739- dmaxp=val;
740- else if((val<0.) && (val<dmaxm))
741- dmaxm=val;
742- }
743- }
744- }
745-
746- /* set all none intersected vertices to max/min distance */
747- for(i=0;i<numvertices;i++){
748- vec_dist_zerolevelset->GetValue(&val,i);
749- if(val==1.) //FIXME: improve check
750- vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
751- else if(val==-1.)
752- vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
753- }
754-
755- /*Assemble vector and serialize */
756- vec_dist_zerolevelset->Assemble();
757- IssmDouble* dist_zerolevelset=vec_dist_zerolevelset->ToMPISerial();
758- InputUpdateFromVectorx(femmodel,dist_zerolevelset,MaskIceLevelsetEnum,VertexSIdEnum);
759-
760- /*Clean up and return*/
761- delete vec_dist_zerolevelset;
762- delete dist_zerolevelset;
763-}/*}}}*/
764 void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
765
766 if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
767Index: ../trunk-jpl/src/c/analyses/analyses.h
768===================================================================
769--- ../trunk-jpl/src/c/analyses/analyses.h (revision 21930)
770+++ ../trunk-jpl/src/c/analyses/analyses.h (revision 21931)
771@@ -33,7 +33,6 @@
772 #include "./SmbAnalysis.h"
773 #include "./SealevelriseAnalysis.h"
774 #include "./MeltingAnalysis.h"
775-#include "./MeshdeformationAnalysis.h"
776 #include "./SmoothAnalysis.h"
777 #include "./StressbalanceAnalysis.h"
778 #include "./StressbalanceSIAAnalysis.h"
779Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.h
780===================================================================
781--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 21930)
782+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 21931)
783@@ -31,7 +31,6 @@
784 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
785 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
786 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
787- void SetDistanceOnIntersectedElements(FemModel* femmodel);
788 void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element);
789 void UpdateConstraints(FemModel* femmodel);
790 };
791Index: ../trunk-jpl/src/c/analyses/EnumToAnalysis.cpp
792===================================================================
793--- ../trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 21930)
794+++ ../trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 21931)
795@@ -109,9 +109,6 @@
796 #ifdef _HAVE_ESA_
797 case EsaAnalysisEnum : return new EsaAnalysis();
798 #endif
799- #ifdef _HAVE_MESHDEFORMATION_
800- case MeshdeformationAnalysisEnum : return new MeshdeformationAnalysis();
801- #endif
802 #ifdef _HAVE_LEVELSET_
803 case LevelsetAnalysisEnum : return new LevelsetAnalysis();
804 #endif
Note: See TracBrowser for help on using the repository browser.