source: issm/oecreview/Archive/17984-18295/ISSM-18223-18224.diff@ 18296

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

Added 17984-18295
File size: 3.9 KB
RevLine 
[18296]1Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18223)
4+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18224)
5@@ -47,8 +47,11 @@
6
7 /*Intermediaries */
8 int num_responses,i;
9- IssmDouble hobs,hu2,weight,NUMx,NUMy,DEN,Jdet;
10- IssmDouble vx,vy,vbar2,nux,nuy,phi,dphi[2];
11+ IssmDouble hobs,hu2,weight,Jdet;
12+ IssmDouble NUMxH2,NUMyH2,DENH2;
13+ IssmDouble NUMxUbar,NUMyUbar,DENUbar;
14+ IssmDouble vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs;
15+ IssmDouble nux,nuy,phi,dphi[2];
16 int *responses = NULL;
17 IssmDouble *xyz_list = NULL;
18
19@@ -67,8 +70,8 @@
20 Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input);
21 Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
22 Input* potential_input = element->GetInput(PotentialEnum); _assert_(potential_input);
23- Input* vx_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input);
24- Input* vy_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input);
25+ Input* vxobs_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
26+ Input* vyobs_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
27 Input* nux_input = element->GetInput(BalancethicknessNuxEnum); _assert_(nux_input);
28 Input* nuy_input = element->GetInput(BalancethicknessNuyEnum); _assert_(nuy_input);
29
30@@ -81,29 +84,43 @@
31 element->NodalFunctions(basis,gauss);
32 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
33
34- vx_input->GetInputValue(&vx,gauss);
35- vy_input->GetInputValue(&vy,gauss);
36+ vxobs_input->GetInputValue(&vxobs,gauss);
37+ vyobs_input->GetInputValue(&vyobs,gauss);
38 nux_input->GetInputValue(&nux,gauss);
39 nuy_input->GetInputValue(&nuy,gauss);
40 potential_input->GetInputValue(&phi,gauss);
41 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
42 thicknessobs_input->GetInputValue(&hobs,gauss);
43
44- vbar2 = (nux*nux*vx*vx + nuy*nuy*vy*vy);
45- hu2 = hobs*hobs*vbar2;
46+ vxobsbar = nux*vxobs;
47+ vyobsbar = nuy*vyobs;
48
49- NUMx = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
50- NUMy = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
51- DEN = vbar2*vbar2+1.e-14;
52+ vbarobs2 = (nux*nux*vxobs*vxobs + nuy*nuy*vyobs*vyobs);
53+ vbarobs = sqrt(vbarobs2);
54+ hu2 = hobs*hobs*vbarobs2;
55
56+ /*H^2 - Hobs^2*/
57+ NUMxH2 = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
58+ NUMyH2 = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
59+ DENH2 = vbarobs2*vbarobs2+1.e-14;
60+
61+ /*Ubar-Ubar_obs*/
62+ NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0];
63+ NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1];
64+ DENUbar = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
65+
66 /*Loop over all requested responses*/
67 for(int resp=0;resp<num_responses;resp++){
68 weights_input->GetInputValue(&weight,gauss,responses[resp]);
69
70 switch(responses[resp]){
71 case Balancethickness2MisfitEnum:
72- //for(i=0;i<numnodes;i++) pe->values[i]+=(NUMx+NUMy)/DEN *weight*Jdet*gauss->weight;
73- for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight;
74+ /*J = (H^2 - Hobs^2)^2*/
75+ //for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight; NOT WORKING
76+ /*J = phi^2*/
77+ //for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; OK
78+ /*J = (ubar - nux*uobs)^2*/
79+ for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
80 break;
81 default:
82 _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
Note: See TracBrowser for help on using the repository browser.