Changeset 18056

Timestamp:

05/25/14 17:19:28 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: changing cost function for bal H 2

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Tria.cpp (modified) (1 diff)
• FemModel.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3176 +3176 @@
 
         /*Intermediaries*/
-        int    vertexpidlist[NUMVERTICES];
-        IssmDouble lambda[NUMVERTICES];
-        IssmDouble gradient_g[NUMVERTICES];
-
-        /*Compute Gradient*/
+        IssmDouble lambda,potential,weight,Jdet;
+        IssmDouble grade_g[NUMVERTICES];
+        IssmDouble basis[NUMVERTICES];
+        IssmDouble xyz_list[NUMVERTICES][3];
+        int        vertexpidlist[NUMVERTICES];
+
+        /*Recover inputs*/
+        Input* lambda_input  = inputs->GetInput(AdjointEnum);                            _assert_(lambda_input);
+        Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+
+        /* Get node coordinates and dof list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         GradientIndexing(&vertexpidlist[0],control_index);
-        GetInputListOnVertices(&lambda[0],AdjointEnum);
-        for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=lambda[i];
-
-        gradient->SetValues(NUMVERTICES,vertexpidlist,gradient_g,INS_VAL);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=new GaussTria(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctions(&basis[0],gauss);
+                weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
+                lambda_input->GetInputValue(&lambda,gauss);
+
+                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+                for(int i=0;i<NUMVERTICES;i++){
+                        grade_g[i]+= - weight*Jdet*gauss->weight*basis[i]*lambda;
+                }
+        }
+        gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
+
+        /*Clean up and return*/
+        delete gauss;
 }
 /*}}}*/

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

@@ -1208 +1208 @@
         IssmDouble J_sum;
 
-        IssmDouble  weight,thicknessobs,thickness;
+        IssmDouble  weight,thicknessobs,thickness,potential;
         IssmDouble  Jdet;
         IssmDouble* xyz_list = NULL;
@@ -1224 +1224 @@
                 Input* thickness_input   =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
                 Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
+                Input* potential_input   =element->GetInput(PotentialEnum);                          _assert_(potential_input);
 
                 /* Start  looping on the number of gaussian points: */
@@ -1238 +1239 @@
                         thickness_input->GetInputValue(&thickness,gauss);
                         thicknessobs_input->GetInputValue(&thicknessobs,gauss);
+                        potential_input->GetInputValue(&potential,gauss);
                         gauss->GaussPoint(ig);
 
-                        J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
+                        //J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
+                        J +=.5*potential*potential*weight*Jdet*gauss->weight;
                 }