Changeset 18222

Timestamp:

07/07/14 20:23:23 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on other cost functions

File:

• issm/trunk-jpl/src/c/classes/FemModel.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1207 +1207 @@
 
         IssmDouble  weight,thicknessobs,thickness,potential;
+        IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nux,nuy;
         IssmDouble  Jdet;
         IssmDouble* xyz_list = NULL;
@@ -1223 +1224 @@
                 Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
                 Input* potential_input   =element->GetInput(PotentialEnum);                          _assert_(potential_input);
+                Input* vxobs_input       =element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
+                Input* vyobs_input       =element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
+                Input* vx_input          =element->GetInput(VxEnum); _assert_(vx_input);
+                Input* vy_input          =element->GetInput(VyEnum); _assert_(vy_input);
+                Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
+                Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
 
                 /* Start  looping on the number of gaussian points: */
@@ -1238 +1245 @@
                         thicknessobs_input->GetInputValue(&thicknessobs,gauss);
                         potential_input->GetInputValue(&potential,gauss);
-                        gauss->GaussPoint(ig);
-
+                        vxobs_input->GetInputValue(&vxobs,gauss);
+                        vyobs_input->GetInputValue(&vyobs,gauss);
+                        vx_input->GetInputValue(&vxbar,gauss);
+                        vy_input->GetInputValue(&vybar,gauss);
+                        nux_input->GetInputValue(&nux,gauss);
+                        nuy_input->GetInputValue(&nuy,gauss);
+
+                        vxbarobs = nux*vxobs;
+                        vybarobs = nuy*vyobs;
+
+                        /*J = (H^2 - Hobs^2)^2*/
                         //J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
-                        J +=.5*potential*potential*weight*Jdet*gauss->weight;
+                        /*J = phi^2*/
+                        //J +=.5*potential*potential*weight*Jdet*gauss->weight; OK
+                        /*J = (ubar - nux*uobs)^2*/
+                        J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
                 }