/*!\file Cfdragcoeffabsgrad.cpp
 * \brief: Cfdragcoeffabsgrad Object
 */

/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
   #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "./classes.h"
#include "./ExternalResults/ExternalResult.h"
#include "./ExternalResults/Results.h"
#include "../datastructures/datastructures.h"
#include "./Elements/Element.h"
#include "./Elements/Elements.h"
#include "./FemModel.h"
#include "../modules/SurfaceAreax/SurfaceAreax.h"
#include "../classes/Params/Parameters.h"
#include "../classes/gauss/Gauss.h"
#include "./Inputs/DatasetInput.h"
/*}}}*/

/*Cfdragcoeffabsgrad constructors, destructors :*/
Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(){/*{{{*/

	this->definitionenum = -1;
	this->name = NULL;
	this->weights_enum = UNDEF;
	this->misfit=0;
	this->timepassedflag = false;
}
/*}}}*/
Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(char* in_name, int in_definitionenum, int in_weights_enum, bool in_timepassedflag){/*{{{*/

	this->definitionenum=in_definitionenum;

	this->name		= xNew<char>(strlen(in_name)+1);
	xMemCpy<char>(this->name,in_name,strlen(in_name)+1);

	this->weights_enum=in_weights_enum;
	this->timepassedflag=in_timepassedflag;

	this->misfit=0;
}
/*}}}*/
Cfdragcoeffabsgrad::~Cfdragcoeffabsgrad(){/*{{{*/
	if(this->name)xDelete(this->name);
	this->misfit=0;
}
/*}}}*/
/*Object virtual function resolutoin: */
Object* Cfdragcoeffabsgrad::copy() {/*{{{*/
	Cfdragcoeffabsgrad* mf = new Cfdragcoeffabsgrad(this->name,this->definitionenum, this->weights_enum,this->timepassedflag);
	mf->misfit=this->misfit;
	return (Object*) mf;
}
/*}}}*/
void Cfdragcoeffabsgrad::DeepEcho(void){/*{{{*/
	this->Echo();
}
/*}}}*/
void Cfdragcoeffabsgrad::Echo(void){/*{{{*/
	_printf_(" Cfdragcoeffabsgrad: " << name << " " << this->definitionenum << "\n");
	_printf_("    weights_enum: " << weights_enum << " " << EnumToStringx(weights_enum) << "\n");
	_printf_("	  timepassedflag: "<<timepassedflag<<"\n");
}
/*}}}*/
int Cfdragcoeffabsgrad::Id(void){/*{{{*/
	return -1;
}
/*}}}*/
void Cfdragcoeffabsgrad::Marshall(MarshallHandle* marshallhandle){/*{{{*/
	_error_("not implemented yet!"); 
} 
/*}}}*/
int Cfdragcoeffabsgrad::ObjectEnum(void){/*{{{*/
	return CfdragcoeffabsgradEnum;
}
/*}}}*/
/*Definition virtual function resolutoin: */
int Cfdragcoeffabsgrad::DefinitionEnum(){/*{{{*/
	return this->definitionenum;
}
/*}}}*/
char* Cfdragcoeffabsgrad::Name(){/*{{{*/
	char* name2=xNew<char>(strlen(this->name)+1);
	xMemCpy(name2,this->name,strlen(this->name)+1);

	return name2;
}
/*}}}*/
IssmDouble Cfdragcoeffabsgrad::Response(FemModel* femmodel){/*{{{*/

	/*recover parameters: */
	IssmDouble J=0.;
	IssmDouble J_sum=0.;

	for(Object* & object : femmodel->elements->objects){
		Element* element=xDynamicCast<Element*>(object);
		J+=this->Cfdragcoeffabsgrad_Calculation(element,weights_enum);
	}

	ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
	ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
	J=J_sum;

	timepassedflag = true;
	return J;
}/*}}}*/
IssmDouble Cfdragcoeffabsgrad::Cfdragcoeffabsgrad_Calculation(Element* element, int weights_enum){/*{{{*/

	int        domaintype,numcomponents;
	IssmDouble Jelem=0.;
	IssmDouble misfit,Jdet;
	IssmDouble dp[2],weight;
	IssmDouble* xyz_list = NULL;

	/*Get basal element*/
	if(!element->IsOnBase()) return 0.;

	/*If on water, return 0: */
	if(!element->IsIceInElement()) return 0.;

	/*Get problem dimension*/
	element->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DverticalEnum:   numcomponents   = 1; break;
		case Domain3DEnum:           numcomponents   = 2; break;
		case Domain2DhorizontalEnum: numcomponents   = 2; break;
		default: _error_("not supported yet");
	}

	/*Spawn surface element*/
	Element* basalelement = element->SpawnBasalElement();

	/* Get node coordinates*/
	basalelement->GetVerticesCoordinates(&xyz_list);

	/*Get input if it already exists*/
	DatasetInput *datasetinput = basalelement->GetDatasetInput(definitionenum);  _assert_(datasetinput);
	Input        *drag_input   = basalelement->GetInput(FrictionCoefficientEnum); _assert_(drag_input);

	/* Start  looping on the number of gaussian points: */
	Gauss* gauss=basalelement->NewGauss(2);
	while(gauss->next()){

		/* Get Jacobian determinant: */
		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);

		/*Get all parameters at gaussian point*/
		datasetinput->GetInputValue(&weight,gauss,WeightsSurfaceObservationEnum);
		drag_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);

		/*Add to cost function*/
		Jelem+=weight*.5*dp[0]*dp[0]*Jdet*gauss->weight;
		if(numcomponents==2) Jelem+=weight*.5*dp[1]*dp[1]*Jdet*gauss->weight;
	}

	/*clean up and Return: */
	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
	xDelete<IssmDouble>(xyz_list);
	delete gauss;
	return Jelem;
}/*}}}*/