/*!\file CreateJacobianMatrixx
 * \brief retrieve vector from inputs in elements
 */

#include "./CreateJacobianMatrixx.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../AllocateSystemMatricesx/AllocateSystemMatricesx.h"

void CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,FemModel* femmodel,IssmDouble kmax){

	int      i,connectivity;
	int      numberofdofspernode;
	int      fsize,configuration_type;
	Element *element = NULL;
	Load    *load    = NULL;
	Matrix<IssmDouble>* Jff = NULL;

	/*Checks*/
	_assert_(femmodel && femmodel->nodes && femmodel->elements);

	/*Recover some parameters*/
	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
	femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
	fsize=femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
	numberofdofspernode=femmodel->nodes->MaxNumDofs(configuration_type,GsetEnum);

	/*Initialize Jacobian Matrix*/
	AllocateSystemMatricesx(&Jff,NULL,NULL,NULL,femmodel);

	/*Create and assemble matrix*/
	for(i=0;i<femmodel->elements->Size();i++){
		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
		element->CreateJacobianMatrix(Jff);
	}
	for (i=0;i<femmodel->loads->Size();i++){
		load=(Load*)femmodel->loads->GetObjectByOffset(i);
		if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff);
		if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax);
	}
	Jff->Assemble();

	/*Assign output pointer*/
	*pJff=Jff;

}
