Changeset 11313 for issm/trunk-jpl/src/c/solutions/gradient_core.cpp

Timestamp:

02/03/12 07:55:52 (13 years ago)

Author:

Mathieu Morlighem

Message:

Objective function gradient is now ONE unique vector

File:

• issm/trunk-jpl/src/c/solutions/gradient_core.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/solutions/gradient_core.cpp ¶

@@ -14 +14 @@
 
 void gradient_core(FemModel* femmodel,int step,bool orthogonalize){ 
-        
-        /*parameters: */
-        bool    control_steady;
-        int     num_controls;
-        int    *control_type   = NULL;
-        double *optscal_list   = NULL;
-        double  optscal=1.e+100,norm_grad;
 
         /*Intermediaries*/
-        Vec new_gradient=NULL;
-        Vec gradient=NULL;
-        Vec old_gradient=NULL;
+        double  norm_inf;
+        double *norm_list    = NULL;
+        Vec     new_gradient = NULL;
+        Vec     gradient     = NULL;
+        Vec     old_gradient = NULL;
 
-        /*retrieve parameters:*/
-        femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
-        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-        femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
-        femmodel->parameters->FindParam(&optscal_list,NULL,NULL,InversionGradientScalingEnum);
+        /*Compute gradient*/
+        _printf_(VerboseControl(),"   compute cost function gradient\n");
+        Gradjx(&gradient,&norm_list,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);
 
-        /*Compute and norm gradient of all controls*/
-        for (int i=0;i<num_controls;i++){
-
-                _printf_(VerboseControl(),"   compute gradient of J with respect to %s\n",EnumToStringx(control_type[i]));
-                Gradjx(&gradient,&norm_grad,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,i);
-
-                if (orthogonalize){
-                        _printf_(VerboseControl(),"   orthogonalization\n");
-                        ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]);
-                        Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient);
-                }
-                else{ 
-                        new_gradient=gradient;
-                }
-
-                /*Get scaling factor of current control:*/
-                VecNorm(new_gradient,NORM_INFINITY,&norm_grad);
-                if(norm_grad<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm of J with respect to %s is zero",EnumToStringx(control_type[i]));
-                if(isnan(norm_grad))_error_("||∂J/∂α||∞ = NaN  gradient norm of J with respect to %s is NaN" ,EnumToStringx(control_type[i]));
-                optscal=min(optscal,optscal_list[num_controls*step+i]/norm_grad);
-
-                /*plug back into inputs: */
-                ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],new_gradient);
-                VecFree(&new_gradient);
+        if (orthogonalize){
+                _printf_(VerboseControl(),"   orthogonalization\n");
+                ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);
+                Orthx(&new_gradient,gradient,old_gradient); VecFree(&old_gradient); VecFree(&gradient);
+        }
+        else{ 
+                new_gradient=gradient;
         }
 
+        /*Check that gradient is clean*/
+        VecNorm(new_gradient,NORM_INFINITY,&norm_inf);
+        if(norm_inf<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm is zero");
+        if(isnan(norm_inf))_error_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
+
+        /*plug back into inputs: */
+        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,new_gradient);
+        VecFree(&new_gradient);
+
         /*Scale Gradients*/
-        for (int i=0;i<num_controls;i++) ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i],optscal);
+        ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,norm_list,step);
 
         /*Clean up and return*/
-        xfree((void**)&control_type);
-        xfree((void**)&optscal_list);
+        xfree((void**)&norm_list);
 }