source: issm/oecreview/Archive/18296-19100/ISSM-18568-18569.diff

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

@@ -136 +136 @@
 
         /*Deal with second term*/
         switch(control_type){
-                case BalancethicknessOmegaEnum:  GradientJOmega(element,gradient,control_index); break;
+                case BalancethicknessOmegaEnum:           GradientJOmega(element,gradient,control_index); break;
+                case BalancethicknessThickeningRateEnum:  GradientJdHdt( element,gradient,control_index); break;
                 default: _error_("control type not supported yet: " << EnumToStringx(control_type));
         }
 
@@ -191 +192 @@
         xDelete<int>(vertexpidlist);
         delete gauss;
 }/*}}}*/
+void AdjointBalancethickness2Analysis::GradientJdHdt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Fetch number of vertices for this finite element*/
+        int numvertices = element->GetNumberOfVertices();
+
+        /*Initialize some vectors*/
+        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
+        IssmDouble* lambda        = xNew<IssmDouble>(numvertices);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /*Retrieve all inputs we will be needing: */
+        element->GradientIndexing(&vertexpidlist[0],control_index);
+        element->GetInputListOnVertices(lambda,AdjointEnum);
+        for(int i=0;i<numvertices;i++){
+                ge[i]= - lambda[i];
+                _assert_(!xIsNan<IssmDouble>(ge[i]));
+        }
+        gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(ge);
+        xDelete<IssmDouble>(lambda);
+        xDelete<int>(vertexpidlist);
+}/*}}}*/
 void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
         element->InputUpdateFromSolutionOneDof(solution,AdjointEnum);
 }/*}}}*/

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

@@ -28 +28 @@
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
                 void GradientJOmega(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void GradientJdHdt(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
                 void UpdateConstraints(FemModel* femmodel);
 };