/*!\file CreateJacobianMatrixx
 * \brief: create system matrices (stiffness matrix, loads vector)
 */

#include "./CreateJacobianMatrixx.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../io/io.h"
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"

void CreateJacobianMatrixx(Matrix** pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,IssmDouble kmax){
	
	int      i,connectivity;
	int      numberofdofspernode;
	int      fsize,configuration_type;
	Element *element = NULL;
	Load    *load    = NULL;
	Matrix*  Jff     = NULL;

	/*Checks*/
	_assert_(nodes && elements);

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

	/*Initialize Jacobian Matrix*/
	Jff=new Matrix(fsize,fsize,connectivity,numberofdofspernode);
	
	/*Create and assemble matrix*/
	for(i=0;i<elements->Size();i++){
		element=(Element*)elements->GetObjectByOffset(i);
		element->CreateJacobianMatrix(Jff);
	}
	for (i=0;i<loads->Size();i++){
		load=(Load*)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;
}
