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

../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp

@@ -47 +47 @@
 
         /*Intermediaries */
         int         num_responses,i;
-        IssmDouble  hobs,hu2,weight,NUMx,NUMy,DEN,Jdet;
-        IssmDouble  vx,vy,vbar2,nux,nuy,phi,dphi[2];
+        IssmDouble  hobs,hu2,weight,Jdet;
+        IssmDouble  NUMxH2,NUMyH2,DENH2;
+        IssmDouble  NUMxUbar,NUMyUbar,DENUbar;
+        IssmDouble  vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs;
+        IssmDouble  nux,nuy,phi,dphi[2];
         int        *responses = NULL;
         IssmDouble *xyz_list  = NULL;
 
@@ -67 +70 @@
         Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
         Input* weights_input      = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
         Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
-        Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input);
-        Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input);
+        Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
+        Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
         Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
         Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
 
@@ -81 +84 @@
                 element->NodalFunctions(basis,gauss);
                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
-                vx_input->GetInputValue(&vx,gauss);
-                vy_input->GetInputValue(&vy,gauss);
+                vxobs_input->GetInputValue(&vxobs,gauss);
+                vyobs_input->GetInputValue(&vyobs,gauss);
                 nux_input->GetInputValue(&nux,gauss);
                 nuy_input->GetInputValue(&nuy,gauss);
                 potential_input->GetInputValue(&phi,gauss);
                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
                 thicknessobs_input->GetInputValue(&hobs,gauss);
 
-                vbar2 = (nux*nux*vx*vx + nuy*nuy*vy*vy);
-                hu2 = hobs*hobs*vbar2;
+                vxobsbar = nux*vxobs;
+                vyobsbar = nuy*vyobs;
 
-                NUMx = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
-                NUMy = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
-                DEN = vbar2*vbar2+1.e-14;
+                vbarobs2 = (nux*nux*vxobs*vxobs + nuy*nuy*vyobs*vyobs);
+                vbarobs  = sqrt(vbarobs2);
+                hu2 = hobs*hobs*vbarobs2;
 
+                /*H^2 - Hobs^2*/
+                NUMxH2 = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
+                NUMyH2 = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
+                DENH2 = vbarobs2*vbarobs2+1.e-14;
+
+                /*Ubar-Ubar_obs*/
+                NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0];
+                NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1];
+                DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
+
                 /*Loop over all requested responses*/
                 for(int resp=0;resp<num_responses;resp++){
                         weights_input->GetInputValue(&weight,gauss,responses[resp]);
 
                         switch(responses[resp]){
                                 case Balancethickness2MisfitEnum:
-                                        //for(i=0;i<numnodes;i++) pe->values[i]+=(NUMx+NUMy)/DEN *weight*Jdet*gauss->weight;
-                                        for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight;
+                                        /*J = (H^2 - Hobs^2)^2*/
+                                        //for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight; NOT WORKING
+                                        /*J = phi^2*/
+                                        //for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; OK
+                                        /*J = (ubar - nux*uobs)^2*/
+                                        for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
                                         break;
                                 default:
                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");