//
// Plane4.h
// plane stress (default implemented) 18.7.2002
// plane strain XXX
// axisymmetric XXX
// 4 node element
//

//Gauss points
extern	Gauss22 *DGauss22;


class Plane4
{
private:
	long number;			
	long NbLink[4];			// 4 connections
	long Mat;				// Material
	long R;					// Real constant for thickness
	double Thickness;		//Thickness of the element
	double **StiffMat;		// 
public:
	Plane4(long num, long mat, long real, long *connectivity, Node *DNodes);
	void CalculateMatrix(Node *DNodes, Material **DMat);
	void AssembleMatrix(double **Global,int MaxDOF,Node *DNodes);
	void WriteStiffMatrix(FILE *file);
};

///////////////////////////////////////

void Plane4::CalculateMatrix(Node *DNodes, Material **DMat)
{
	double Jacobian [2][2],JacDet;
	double EMatrixConst;
	double NShFu[4];		// shape function
	double NDerR[4];		// derivative relative to R
	double NDerS[4];		// derivative relative to S

	double **EMatrix=new double*[3];
	double **BMatrix=new double*[3];		// B Matrix
	double **BtMatrix=new double*[8];		// Transpose of B Matrix
	double **TempMatrix=new double*[3];
	double temp02[8][8];

	int i,j,k;


	for (i=0;i<8;i++)
		for(j=0;j<8;j++)
			temp02[i][j]=0.0;				// initialize matrix
		
	for (i=0;i<3;i++)			// 3*8
	{
		BMatrix[i] = new double[8];
		TempMatrix[i] = new double[8];
	}
	
	for (i=0;i<3;i++)			// 3*3
		EMatrix[i] = new double[3];
	
	for (i=0;i<8;i++)			// 8*3
		BtMatrix[i] = new double[3];

	StiffMat= new double*[8];
	for (i=0;i<8;i++)					// 8*8
		StiffMat[i] = new double[8];

	//Calculation of EMatrix
	EMatrixConst=DMat[Mat]->Ex/(1-pow(DMat[Mat]->Nuxy,2));
	
	for (i=0;i<3;i++)
		for (j=0;j<3;j++)
			EMatrix[i][j]=0;

	EMatrix[0][0]=EMatrixConst;
	EMatrix[0][1]=(DMat[Mat]->Nuxy*EMatrixConst);
	EMatrix[1][0]=(DMat[Mat]->Nuxy*EMatrixConst);
	EMatrix[1][1]=EMatrixConst;
	EMatrix[2][2]=(0.5*(1-DMat[Mat]->Nuxy)*EMatrixConst);


	for (i=0;i<4;i++)	{ // where 2 is the order of the integration (2*2 int points)
		// Definition of Ni, and their derivatives
		NShFu[0]= 0.25*(1-DGauss22->Rg[i])*(1-DGauss22->Sg[i]);
		NShFu[1]= 0.25*(1+DGauss22->Rg[i])*(1-DGauss22->Sg[i]);
		NShFu[2]= 0.25*(1+DGauss22->Rg[i])*(1+DGauss22->Sg[i]);
		NShFu[3]= 0.25*(1-DGauss22->Rg[i])*(1+DGauss22->Sg[i]);
			
		NDerR[0]= -0.25*(1-DGauss22->Sg[i]);
		NDerR[1]=  0.25*(1-DGauss22->Sg[i]);
		NDerR[2]=  0.25*(1+DGauss22->Sg[i]);
		NDerR[3]= -0.25*(1+DGauss22->Sg[i]);
				
		NDerS[0]= -0.25*(1-DGauss22->Rg[i]);
		NDerS[1]= -0.25*(1+DGauss22->Rg[i]);
		NDerS[2]=  0.25*(1+DGauss22->Rg[i]);
		NDerS[3]=  0.25*(1-DGauss22->Rg[i]);

		// initialisation of the Jacobian Matrix
		for (j=0;j<2;j++)
			for (k=0;k<2;k++)
				Jacobian[j][k]=0.0;
		
		for(j=0;j<4;j++)
		{
			Jacobian[0][0]+=NDerR[j]*DNodes[(NbLink[j]-1)].nx[0];
			Jacobian[0][1]+=NDerR[j]*DNodes[(NbLink[j]-1)].nx[1];
			Jacobian[1][0]+=NDerS[j]*DNodes[(NbLink[j]-1)].nx[0]; 
			Jacobian[1][1]+=NDerS[j]*DNodes[(NbLink[j]-1)].nx[1];
		}
		
		JacDet=Jacobian[0][0]*Jacobian[1][1]-Jacobian[1][0]*Jacobian[0][1];
		if(JacDet==0.0)
			nerror("Element Jacobian must be greater than 0.");

		// Calculation of BMatrix[][]
		for(j=0;j<3;j++)
			for(k=0;k<8;k++)
				BMatrix[j][k] = 0.0;
		

		for (j=0;j<4;j++)
		{
			BMatrix[0][j*2]   = Jacobian[1][1]*NDerR[j]-Jacobian[0][1]*NDerS[j];
			BMatrix[1][j*2+1] = -Jacobian[1][0]*NDerR[j]+Jacobian[0][0]*NDerS[j];
			BMatrix[2][j*2]   = BMatrix[1][j*2+1];
			BMatrix[2][j*2+1] = BMatrix[0][j*2];
		}

		for (j=0;j<3;j++)
			for(k=0;k<8;k++)
				BMatrix[j][k] = BMatrix[j][k]/JacDet;
		
		// Calculation of the Transpose of BMatrix[][]
		matrix_transpose(BtMatrix,BMatrix,8,3);
		
		matrix_multiply(EMatrix, BMatrix, TempMatrix,3,3,8);
		matrix_multiply(BtMatrix, TempMatrix, StiffMat,8,3,8);
		
		// Gauss integration: multiplied by 1 for the weight
		for(j=0;j<8;j++)
			for(k=0;k<8;k++)
				temp02[j][k]+=StiffMat[j][k]*JacDet;
		}
	for(j=0;j<8;j++)
		for(k=0;k<8;k++)
			StiffMat[j][k]=temp02[j][k];
}

///////////////////////////////////////

Plane4::Plane4(long num, long mat, long real, long *connectivity, Node *DNodes)
{
	int i;

	number = num;
	Mat = mat;
	R = real;

	for(i=0;i<4;i++)	NbLink[i] = connectivity[i];

	for(i=0;i<4;i++)	{
		DNodes[(NbLink[i]-1)].ActiveDOF[0]=1;	// x DOF active node i
		DNodes[(NbLink[i]-1)].ActiveDOF[1]=1;	// y DOF active
		}
}

///////////////////////////////////////

void Plane4::AssembleMatrix(double **Global, int MaxDOF,Node *DNodes)
{
	int i,j,k,l;

	for(i=0;i<4;i++)				// step for column (3 Nodes 2x2)
		for(j=0;j<4;j++)			// step for row (3 Nodes 2x2)
			for(k=0;k<2;k++)		// step for line (row) (2 DOF)
				for(l=0;l<2;l++)	// step for line (column) (2 DOF)
					Global[(DNodes[(NbLink[i]-1)].CumSum+k)][(DNodes[NbLink[j]-1].CumSum+l)]+=StiffMat[(i*2)+k][(j*2)+l];
}

///////////////////////////////////////

void Plane4::WriteStiffMatrix(FILE *file)
{
	int i,j;

	fprintf(file,"Element matrix %7d\n",number);
	for (i=0;i<8;i++)
	{
		for (j=0;j<8;j++)	{
			fprintf(file,"%13.3E ",StiffMat[i][j]);
			}
		fprintf(file,"\n");
		}
	fprintf(file,"\n");
}