source: issm/oecreview/Archive/17802-17983/ISSM-17974-17975.diff@ 17986

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

New release
File size: 28.3 KB
RevLine 
[17986]1Index: ../trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 17974)
4+++ ../trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 17975)
5@@ -9,18 +9,14 @@
6
7 void SurfaceLogVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
8
9- /*Intermediary*/
10- int i;
11- Element* element=NULL;
12-
13 /*output: */
14- IssmDouble J=0;
15+ IssmDouble J=0.;
16 IssmDouble J_sum;
17
18 /*Compute Misfit: */
19- for (i=0;i<elements->Size();i++){
20- element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
21- J+=element->SurfaceLogVelMisfit();
22+ for(int i=0;i<elements->Size();i++){
23+ Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
24+ J+=SurfaceLogVelMisfit(element);
25 }
26
27 /*Sum all J from all cpus of the cluster:*/
28@@ -31,3 +27,102 @@
29 /*Assign output pointers: */
30 *pJ=J;
31 }
32+
33+IssmDouble SurfaceLogVelMisfit(Element* element){
34+
35+ int domaintype,numcomponents;
36+ bool surfaceintegral;
37+ IssmDouble Jelem=0.;
38+ IssmDouble epsvel=2.220446049250313e-16;
39+ IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
40+ IssmDouble velocity_mag,obs_velocity_mag;
41+ IssmDouble misfit,Jdet;
42+ IssmDouble vx,vy,vxobs,vyobs,weight;
43+ IssmDouble* xyz_list = NULL;
44+
45+ /*Get basal element*/
46+ if(!element->IsOnSurface()) return 0.;
47+
48+ /*If on water, return 0: */
49+ if(!element->IsIceInElement()) return 0.;
50+
51+ /*Get problem dimension*/
52+ element->FindParam(&domaintype,DomainTypeEnum);
53+ switch(domaintype){
54+ case Domain2DverticalEnum:
55+ surfaceintegral = true;
56+ numcomponents = 1;
57+ break;
58+ case Domain3DEnum:
59+ surfaceintegral = true;
60+ numcomponents = 2;
61+ break;
62+ case Domain2DhorizontalEnum:
63+ surfaceintegral = false;
64+ numcomponents = 2;
65+ break;
66+ default: _error_("not supported yet");
67+ }
68+
69+ /*Spawn surface element*/
70+ Element* topelement = element->SpawnTopElement();
71+
72+ /* Get node coordinates*/
73+ topelement->GetVerticesCoordinates(&xyz_list);
74+
75+ /*Retrieve all inputs we will be needing: */
76+ Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
77+ Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input);
78+ Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
79+ Input* vy_input = NULL;
80+ Input* vyobs_input = NULL;
81+ if(numcomponents==2){
82+ vy_input =topelement->GetInput(VyEnum); _assert_(vy_input);
83+ vyobs_input =topelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
84+ }
85+
86+ /* Start looping on the number of gaussian points: */
87+ Gauss* gauss=topelement->NewGauss(4);
88+ for(int ig=gauss->begin();ig<gauss->end();ig++){
89+
90+ gauss->GaussPoint(ig);
91+
92+ /* Get Jacobian determinant: */
93+ topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
94+
95+ /*Get all parameters at gaussian point*/
96+ weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
97+ vx_input->GetInputValue(&vx,gauss);
98+ vxobs_input->GetInputValue(&vxobs,gauss);
99+ if(numcomponents==2){
100+ vy_input->GetInputValue(&vy,gauss);
101+ vyobs_input->GetInputValue(&vyobs,gauss);
102+ }
103+
104+ /*Compute SurfaceLogVelMisfit:
105+ * [ vel + eps ] 2
106+ * J = 4 \bar{v}^2 | log ( ----------- ) |
107+ * [ vel + eps ]
108+ * obs
109+ */
110+ if(numcomponents==1){
111+ velocity_mag =fabs(vx)+epsvel;
112+ obs_velocity_mag=fabs(vxobs)+epsvel;
113+ }
114+ else{
115+ velocity_mag =sqrt(vx*vx+vy*vy)+epsvel;
116+ obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel;
117+ }
118+
119+ misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
120+
121+ /*Add to cost function*/
122+ Jelem+=misfit*weight*Jdet*gauss->weight;
123+ }
124+
125+ /*clean up and Return: */
126+ if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
127+ xDelete<IssmDouble>(xyz_list);
128+ delete gauss;
129+ return Jelem;
130+}
131Index: ../trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h
132===================================================================
133--- ../trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h (revision 17974)
134+++ ../trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h (revision 17975)
135@@ -9,5 +9,6 @@
136
137 /* local prototypes: */
138 void SurfaceLogVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
139+IssmDouble SurfaceLogVelMisfit(Element* element);
140
141 #endif
142Index: ../trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp
143===================================================================
144--- ../trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 17974)
145+++ ../trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 17975)
146@@ -9,18 +9,14 @@
147
148 void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
149
150- /*Intermediary*/
151- int i;
152- Element* element=NULL;
153-
154 /*output: */
155- IssmDouble J=0;
156+ IssmDouble J=0.;
157 IssmDouble J_sum;
158
159 /*Compute Misfit: */
160- for (i=0;i<elements->Size();i++){
161- element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
162- J+=element->SurfaceRelVelMisfit();
163+ for(int i=0;i<elements->Size();i++){
164+ Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
165+ J+=SurfaceRelVelMisfit(element);
166 }
167
168 /*Sum all J from all cpus of the cluster:*/
169@@ -31,3 +27,102 @@
170 /*Assign output pointers: */
171 *pJ=J;
172 }
173+
174+IssmDouble SurfaceRelVelMisfit(Element* element){
175+
176+ int domaintype,numcomponents;
177+ bool surfaceintegral;
178+ IssmDouble Jelem=0.;
179+ IssmDouble misfit,Jdet,scalex,scaley;
180+ IssmDouble epsvel=2.220446049250313e-16;
181+ IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
182+ IssmDouble vx,vy,vxobs,vyobs,weight;
183+ IssmDouble* xyz_list = NULL;
184+
185+ /*Get basal element*/
186+ if(!element->IsOnSurface()) return 0.;
187+
188+ /*If on water, return 0: */
189+ if(!element->IsIceInElement()) return 0.;
190+
191+ /*Get problem dimension*/
192+ element->FindParam(&domaintype,DomainTypeEnum);
193+ switch(domaintype){
194+ case Domain2DverticalEnum:
195+ surfaceintegral = true;
196+ numcomponents = 1;
197+ break;
198+ case Domain3DEnum:
199+ surfaceintegral = true;
200+ numcomponents = 2;
201+ break;
202+ case Domain2DhorizontalEnum:
203+ surfaceintegral = false;
204+ numcomponents = 2;
205+ break;
206+ default: _error_("not supported yet");
207+ }
208+
209+ /*Spawn surface element*/
210+ Element* topelement = element->SpawnTopElement();
211+
212+ /* Get node coordinates*/
213+ topelement->GetVerticesCoordinates(&xyz_list);
214+
215+ /*Retrieve all inputs we will be needing: */
216+ Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
217+ Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input);
218+ Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
219+ Input* vy_input = NULL;
220+ Input* vyobs_input = NULL;
221+ if(numcomponents==2){
222+ vy_input =topelement->GetInput(VyEnum); _assert_(vy_input);
223+ vyobs_input =topelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
224+ }
225+
226+ /* Start looping on the number of gaussian points: */
227+ Gauss* gauss=topelement->NewGauss(4);
228+ for(int ig=gauss->begin();ig<gauss->end();ig++){
229+
230+ gauss->GaussPoint(ig);
231+
232+ /* Get Jacobian determinant: */
233+ topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
234+
235+ /*Get all parameters at gaussian point*/
236+ weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
237+ vx_input->GetInputValue(&vx,gauss);
238+ vxobs_input->GetInputValue(&vxobs,gauss);
239+ if(numcomponents==2){
240+ vy_input->GetInputValue(&vy,gauss);
241+ vyobs_input->GetInputValue(&vyobs,gauss);
242+ }
243+
244+ /*Compute SurfaceRelVelMisfit:
245+ *
246+ * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
247+ * J = --- | ------------- (u - u ) + ------------- (v - v ) |
248+ * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
249+ * obs obs
250+ */
251+
252+ if(numcomponents==2){
253+ scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
254+ scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
255+ misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
256+ }
257+ else{
258+ scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
259+ misfit=0.5*(scalex*pow((vx-vxobs),2));
260+ }
261+
262+ /*Add to cost function*/
263+ Jelem+=misfit*weight*Jdet*gauss->weight;
264+ }
265+
266+ /*clean up and Return: */
267+ if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
268+ xDelete<IssmDouble>(xyz_list);
269+ delete gauss;
270+ return Jelem;
271+}
272Index: ../trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h
273===================================================================
274--- ../trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h (revision 17974)
275+++ ../trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h (revision 17975)
276@@ -9,5 +9,6 @@
277
278 /* local prototypes: */
279 void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
280+IssmDouble SurfaceRelVelMisfit(Element* element);
281
282 #endif
283Index: ../trunk-jpl/src/c/classes/Elements/Element.h
284===================================================================
285--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17974)
286+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17975)
287@@ -276,9 +276,6 @@
288
289 virtual void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index)=0;
290 virtual IssmDouble ThicknessAbsMisfit(void)=0;
291- virtual IssmDouble SurfaceAbsVelMisfit(void)=0;
292- virtual IssmDouble SurfaceRelVelMisfit(void)=0;
293- virtual IssmDouble SurfaceLogVelMisfit(void)=0;
294 virtual IssmDouble SurfaceLogVxVyMisfit(void)=0;
295 virtual IssmDouble SurfaceAverageVelMisfit(void)=0;
296 virtual IssmDouble ThicknessAbsGradient(void)=0;
297@@ -286,7 +283,6 @@
298 virtual IssmDouble ThicknessAcrossGradient(void)=0;
299 virtual IssmDouble BalancethicknessMisfit(void)=0;
300 virtual IssmDouble RheologyBbarAbsGradient(void)=0;
301- virtual IssmDouble DragCoefficientAbsGradient(void)=0;
302 virtual void ControlInputGetGradient(Vector<IssmDouble>* gradient,int enum_type,int control_index)=0;
303 virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
304 virtual void ControlInputScaleGradient(int enum_type, IssmDouble scale)=0;
305Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
306===================================================================
307--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17974)
308+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17975)
309@@ -3570,65 +3570,6 @@
310 return Jelem;
311 }
312 /*}}}*/
313-IssmDouble Tria::SurfaceLogVelMisfit(void){/*{{{*/
314-
315- IssmDouble Jelem=0.;
316- IssmDouble misfit,Jdet;
317- IssmDouble epsvel=2.220446049250313e-16;
318- IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
319- IssmDouble velocity_mag,obs_velocity_mag;
320- IssmDouble xyz_list[NUMVERTICES][3];
321- IssmDouble vx,vy,vxobs,vyobs,weight;
322- GaussTria *gauss=NULL;
323-
324- /*If on water, return 0: */
325- if(!IsIceInElement())return 0;
326-
327- /* Get node coordinates and dof list: */
328- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
329-
330- /*Retrieve all inputs we will be needing: */
331- Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
332- Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
333- Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
334- Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
335- Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
336-
337- /* Start looping on the number of gaussian points: */
338- gauss=new GaussTria(4);
339- for(int ig=gauss->begin();ig<gauss->end();ig++){
340-
341- gauss->GaussPoint(ig);
342-
343- /* Get Jacobian determinant: */
344- GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
345-
346- /*Get all parameters at gaussian point*/
347- weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
348- vx_input->GetInputValue(&vx,gauss);
349- vy_input->GetInputValue(&vy,gauss);
350- vxobs_input->GetInputValue(&vxobs,gauss);
351- vyobs_input->GetInputValue(&vyobs,gauss);
352-
353- /*Compute SurfaceLogVelMisfit:
354- * 4 [ vel + eps ] 2
355- * J = 4 \bar{v}^2 | log ( ----------- ) |
356- * [ vel + eps ]
357- * obs
358- */
359- velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel;
360- obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
361- misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
362-
363- /*Add to cost function*/
364- Jelem+=misfit*weight*Jdet*gauss->weight;
365- }
366-
367- /*clean-up and Return: */
368- delete gauss;
369- return Jelem;
370-}
371-/*}}}*/
372 IssmDouble Tria::SurfaceLogVxVyMisfit(void){/*{{{*/
373
374 IssmDouble Jelem=0, S=0;
375@@ -3688,121 +3629,6 @@
376 return Jelem;
377 }
378 /*}}}*/
379-IssmDouble Tria::SurfaceAbsVelMisfit(void){/*{{{*/
380-
381- IssmDouble Jelem=0;
382- IssmDouble misfit,Jdet;
383- IssmDouble vx,vy,vxobs,vyobs,weight;
384- IssmDouble xyz_list[NUMVERTICES][3];
385- GaussTria *gauss=NULL;
386-
387- /*If on water, return 0: */
388- if(!IsIceInElement())return 0;
389-
390- /* Get node coordinates and dof list: */
391- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
392-
393- /*Retrieve all inputs we will be needing: */
394- Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
395- Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
396- Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
397- Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
398- Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
399-
400- /* Start looping on the number of gaussian points: */
401- gauss=new GaussTria(2);
402- for(int ig=gauss->begin();ig<gauss->end();ig++){
403-
404- gauss->GaussPoint(ig);
405-
406- /* Get Jacobian determinant: */
407- GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
408-
409- /*Get all parameters at gaussian point*/
410- weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
411- vx_input->GetInputValue(&vx,gauss);
412- vy_input->GetInputValue(&vy,gauss);
413- vxobs_input->GetInputValue(&vxobs,gauss);
414- vyobs_input->GetInputValue(&vyobs,gauss);
415-
416- /*Compute SurfaceAbsVelMisfitEnum:
417- *
418- * 1 [ 2 2 ]
419- * J = --- | (u - u ) + (v - v ) |
420- * 2 [ obs obs ]
421- *
422- */
423- misfit=0.5*( pow(vx-vxobs,2) + pow(vy-vyobs,2) );
424-
425- /*Add to cost function*/
426- Jelem+=misfit*weight*Jdet*gauss->weight;
427- }
428-
429- /*clean up and Return: */
430- delete gauss;
431- return Jelem;
432-}
433-/*}}}*/
434-IssmDouble Tria::SurfaceRelVelMisfit(void){/*{{{*/
435-
436- IssmDouble Jelem=0;
437- IssmDouble scalex=1,scaley=1;
438- IssmDouble misfit,Jdet;
439- IssmDouble epsvel=2.220446049250313e-16;
440- IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
441- IssmDouble vx,vy,vxobs,vyobs,weight;
442- IssmDouble xyz_list[NUMVERTICES][3];
443- GaussTria *gauss=NULL;
444-
445- /*If on water, return 0: */
446- if(!IsIceInElement())return 0;
447-
448- /* Get node coordinates and dof list: */
449- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
450-
451- /*Retrieve all inputs we will be needing: */
452- Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
453- Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
454- Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
455- Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
456- Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
457-
458- /* Start looping on the number of gaussian points: */
459- gauss=new GaussTria(4);
460- for(int ig=gauss->begin();ig<gauss->end();ig++){
461-
462- gauss->GaussPoint(ig);
463-
464- /* Get Jacobian determinant: */
465- GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
466-
467- /*Get all parameters at gaussian point*/
468- weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
469- vx_input->GetInputValue(&vx,gauss);
470- vy_input->GetInputValue(&vy,gauss);
471- vxobs_input->GetInputValue(&vxobs,gauss);
472- vyobs_input->GetInputValue(&vyobs,gauss);
473-
474- /*Compute SurfaceRelVelMisfit:
475- *
476- * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
477- * J = --- | ------------- (u - u ) + ------------- (v - v ) |
478- * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
479- * obs obs
480- */
481- scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
482- scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
483- misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
484-
485- /*Add to cost function*/
486- Jelem+=misfit*weight*Jdet*gauss->weight;
487- }
488-
489- /*clean up and Return: */
490- delete gauss;
491- return Jelem;
492-}
493-/*}}}*/
494 IssmDouble Tria::ThicknessAbsGradient(void){/*{{{*/
495
496 /* Intermediaries */
497@@ -3988,48 +3814,6 @@
498 return Jelem;
499 }
500 /*}}}*/
501-IssmDouble Tria::DragCoefficientAbsGradient(void){/*{{{*/
502-
503- /* Intermediaries */
504- IssmDouble Jelem = 0;
505- IssmDouble weight;
506- IssmDouble Jdet;
507- IssmDouble xyz_list[NUMVERTICES][3];
508- IssmDouble dp[NDOF2];
509- GaussTria *gauss = NULL;
510-
511- /*retrieve parameters and inputs*/
512-
513- /*If on water, return 0: */
514- if(!IsIceInElement()) return 0;
515-
516- /*Retrieve all inputs we will be needing: */
517- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
518- Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
519- Input* drag_input =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input);
520-
521- /* Start looping on the number of gaussian points: */
522- gauss=new GaussTria(2);
523- for(int ig=gauss->begin();ig<gauss->end();ig++){
524-
525- gauss->GaussPoint(ig);
526-
527- /* Get Jacobian determinant: */
528- GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
529-
530- /*Get all parameters at gaussian point*/
531- weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);
532- drag_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
533-
534- /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
535- Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
536- }
537-
538- /*Clean up and return*/
539- delete gauss;
540- return Jelem;
541-}
542-/*}}}*/
543 void Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
544
545 int vertexpidlist[NUMVERTICES];
546Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
547===================================================================
548--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17974)
549+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17975)
550@@ -136,7 +136,6 @@
551 void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
552 #endif
553
554- IssmDouble DragCoefficientAbsGradient(void);
555 void GradientIndexing(int* indexing,int control_index);
556 void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
557 void GradjBGradient(Vector<IssmDouble>* gradient,int control_index);
558@@ -158,13 +157,10 @@
559 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
560 IssmDouble RheologyBbarAbsGradient(void);
561 IssmDouble ThicknessAbsMisfit(void);
562- IssmDouble SurfaceAbsVelMisfit(void);
563 IssmDouble ThicknessAbsGradient(void);
564 IssmDouble ThicknessAlongGradient(void);
565 IssmDouble ThicknessAcrossGradient(void);
566 IssmDouble BalancethicknessMisfit(void);
567- IssmDouble SurfaceRelVelMisfit(void);
568- IssmDouble SurfaceLogVelMisfit(void);
569 IssmDouble SurfaceLogVxVyMisfit(void);
570 IssmDouble SurfaceAverageVelMisfit(void);
571 void InputControlUpdate(IssmDouble scalar,bool save_parameter);
572Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
573===================================================================
574--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17974)
575+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17975)
576@@ -3364,78 +3364,6 @@
577 }
578 }
579 /*}}}*/
580-IssmDouble Penta::SurfaceAbsVelMisfit(void){/*{{{*/
581-
582- int approximation;
583- IssmDouble J;
584- Tria* tria=NULL;
585-
586- /*retrieve inputs :*/
587- inputs->GetInputValue(&approximation,ApproximationEnum);
588-
589- /*If on water, return 0: */
590- if(!IsIceInElement())return 0;
591-
592- /*Bail out if this element if:
593- * -> Non SSA and not on the surface
594- * -> SSA (2d model) and not on bed) */
595- if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
596- return 0;
597- }
598- else if (approximation==SSAApproximationEnum){
599-
600- /*This element should be collapsed into a tria element at its base. Create this tria element,
601- * and compute SurfaceAbsVelMisfit*/
602- tria=(Tria*)SpawnTria(0,1,2);
603- J=tria->SurfaceAbsVelMisfit();
604- delete tria->material; delete tria;
605- return J;
606- }
607- else{
608-
609- tria=(Tria*)SpawnTria(3,4,5);
610- J=tria->SurfaceAbsVelMisfit();
611- delete tria->material; delete tria;
612- return J;
613- }
614-}
615-/*}}}*/
616-IssmDouble Penta::SurfaceLogVelMisfit(void){/*{{{*/
617-
618- int approximation;
619- IssmDouble J;
620- Tria* tria=NULL;
621-
622- /*retrieve inputs :*/
623- inputs->GetInputValue(&approximation,ApproximationEnum);
624-
625- /*If on water, return 0: */
626- if(!IsIceInElement())return 0;
627-
628- /*Bail out if this element if:
629- * -> Non SSA and not on the surface
630- * -> SSA (2d model) and not on bed) */
631- if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
632- return 0;
633- }
634- else if (approximation==SSAApproximationEnum){
635-
636- /*This element should be collapsed into a tria element at its base. Create this tria element,
637- * and compute SurfaceLogVelMisfit*/
638- tria=(Tria*)SpawnTria(0,1,2); //lower face is 0, upper face is 1.
639- J=tria->SurfaceLogVelMisfit();
640- delete tria->material; delete tria;
641- return J;
642- }
643- else{
644-
645- tria=(Tria*)SpawnTria(3,4,5); //lower face is 0, upper face is 1.
646- J=tria->SurfaceLogVelMisfit();
647- delete tria->material; delete tria;
648- return J;
649- }
650-}
651-/*}}}*/
652 IssmDouble Penta::SurfaceLogVxVyMisfit(void){/*{{{*/
653
654 IssmDouble J;
655@@ -3474,42 +3402,6 @@
656 }
657 }
658 /*}}}*/
659-IssmDouble Penta::SurfaceRelVelMisfit(void){/*{{{*/
660-
661- int approximation;
662- IssmDouble J;
663- Tria* tria=NULL;
664-
665- /*retrieve inputs :*/
666- inputs->GetInputValue(&approximation,ApproximationEnum);
667-
668- /*If on water, return 0: */
669- if(!IsIceInElement())return 0;
670-
671- /*Bail out if this element if:
672- * -> Non SSA and not on the surface
673- * -> SSA (2d model) and not on bed) */
674- if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
675- return 0;
676- }
677- else if (approximation==SSAApproximationEnum){
678-
679- /*This element should be collapsed into a tria element at its base. Create this tria element,
680- * and compute SurfaceRelVelMisfit*/
681- tria=(Tria*)SpawnTria(0,1,2);
682- J=tria->SurfaceRelVelMisfit();
683- delete tria->material; delete tria;
684- return J;
685- }
686- else{
687-
688- tria=(Tria*)SpawnTria(3,4,5);
689- J=tria->SurfaceRelVelMisfit();
690- delete tria->material; delete tria;
691- return J;
692- }
693-}
694-/*}}}*/
695 IssmDouble Penta::ThicknessAbsGradient(void){/*{{{*/
696
697 _error_("Not implemented yet");
698@@ -3534,20 +3426,6 @@
699 return J;
700 }
701 /*}}}*/
702-IssmDouble Penta::DragCoefficientAbsGradient(void){/*{{{*/
703-
704- IssmDouble J;
705- Tria* tria=NULL;
706-
707- /*If on water, on shelf or not on bed, skip: */
708- if(!IsIceInElement()|| IsFloating() || !IsOnBase()) return 0;
709-
710- tria=(Tria*)SpawnTria(0,1,2); //lower face is 0, upper face is 1
711- J=tria->DragCoefficientAbsGradient();
712- delete tria->material; delete tria;
713- return J;
714-}
715-/*}}}*/
716 IssmDouble Penta::RheologyBbarAbsGradient(void){/*{{{*/
717
718 IssmDouble J;
719Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
720===================================================================
721--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17974)
722+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17975)
723@@ -131,7 +131,6 @@
724 void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
725 #endif
726
727- IssmDouble DragCoefficientAbsGradient(void);
728 void GradientIndexing(int* indexing,int control_index);
729 void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
730 void GradjDragSSA(Vector<IssmDouble>* gradient,int control_index);
731@@ -148,9 +147,6 @@
732 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
733 IssmDouble RheologyBbarAbsGradient(void);
734 IssmDouble ThicknessAbsMisfit(void);
735- IssmDouble SurfaceAbsVelMisfit(void);
736- IssmDouble SurfaceRelVelMisfit(void);
737- IssmDouble SurfaceLogVelMisfit(void);
738 IssmDouble SurfaceLogVxVyMisfit(void);
739 IssmDouble SurfaceAverageVelMisfit(void);
740 IssmDouble ThicknessAbsGradient(void);
741Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
742===================================================================
743--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17974)
744+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17975)
745@@ -167,7 +167,6 @@
746 void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
747 #endif
748
749- IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");};
750 void GradientIndexing(int* indexing,int control_index);
751 void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){_error_("not implemented yet");};
752 void GradjBGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
753@@ -189,13 +188,10 @@
754 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
755 IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");};
756 IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");};
757- IssmDouble SurfaceAbsVelMisfit(void){_error_("not implemented yet");};
758 IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");};
759 IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");};
760 IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");};
761 IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");};
762- IssmDouble SurfaceRelVelMisfit(void){_error_("not implemented yet");};
763- IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");};
764 IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");};
765 IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");};
766 void InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
767Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
768===================================================================
769--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17974)
770+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17975)
771@@ -193,13 +193,10 @@
772 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
773 IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");};
774 IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");};
775- IssmDouble SurfaceAbsVelMisfit(void){_error_("not implemented yet");};
776 IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");};
777 IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");};
778 IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");};
779 IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");};
780- IssmDouble SurfaceRelVelMisfit(void){_error_("not implemented yet");};
781- IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");};
782 IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");};
783 IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");};
784 void InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
Note: See TracBrowser for help on using the repository browser.