/*!\file Gradjx
 * \brief: compute inverse method gradient
 */

#include "./Gradjx.h"
#include "../modules.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"

void Gradjx(Vec* pgradient,double** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){

	int     i,j,numberofvertices;
	int     num_controls;
	double  norm_inf;
	double *norm_list       = NULL;
	int    *control_type    = NULL;
	Vec     gradient        = NULL;
	Vec    *gradient_list   = NULL;
	
	/*retrieve some parameters: */
	parameters->FindParam(&num_controls,InversionNumControlParametersEnum);   _assert_(num_controls);
	parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
	numberofvertices=vertices->NumberOfVertices();

	/*Allocate gradient_list */
	gradient_list = (Vec*)xmalloc(num_controls*sizeof(Vec));
	norm_list = (double*)xmalloc(num_controls*sizeof(double));
	for(i=0;i<num_controls;i++){
		gradient_list[i]=NewVec(num_controls*numberofvertices);
	}
	gradient=NewVec(num_controls*numberofvertices);

	/*Compute all gradient_list*/
	for(i=0;i<num_controls;i++){

		for(j=0;j<elements->Size();j++){
			Element* element=(Element*)elements->GetObjectByOffset(j);
			element->Gradj(gradient_list[i],control_type[i],i);
		}

		VecAssemblyBegin(gradient_list[i]);
		VecAssemblyEnd(gradient_list[i]);

		VecNorm(gradient_list[i],NORM_INFINITY,&norm_list[i]);
	}

	/*Add all gradient_list together*/
	for(i=0;i<num_controls;i++){
		VecAXPY(gradient,1.,gradient_list[i]);
		VecFree(&gradient_list[i]);
	}

	/*Check that gradient is clean*/
	VecNorm(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");

	/*Clean-up and assign output pointer*/
	if(pnorm_list){
		*pnorm_list=norm_list;
	}
	else{
		xfree((void**)&norm_list);
	}
	if(pgradient)  *pgradient=gradient;
	xfree((void**)&gradient_list);
	xfree((void**)&control_type);
}
