//
// 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 Mat;				// Material
	long R;					// Real constant for thickness
	double Thickness;		//Thickness of the element
	long NbLink[4];					// 4 connections
	double **StiffMat;		// 
	double **BMatrix;				//3*8 Matrix
	double **EMatrix;				//3*3 
	double **BtMatrix;
	double **TempMatrix;
	
public:
	Plane4(long num, long mat, long real, long *connectivity, Node *DNodes);
	void CalculateMatrix(Node *DNodes, Material **DMat, long Loop);
	void AssembleMatrix(double **Global,int MaxDOF,Node *DNodes);
	void WriteStiffMatrix(FILE *file);
	void CalculateStress(Node *DNodes, Gauss22 *DGauss22);
	double **StrainMatrix;		//3 element array
	double **StressMatrix;		//3 element array
	double ***GaussB; 			//3*8 matrices to store BMatrix for each Gauss point 
};

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

Plane4::Plane4(long num, long mat, long real, long *connectivity, Node *DNodes)
{
	int i,j;
	
	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
		DNodes[(NbLink[i]-1)].connected++;
	}
	
	BtMatrix=new double*[8];		// Transpose of B Matrix
	TempMatrix=new double*[3];
	
	for (i=0;i<3;i++)			// 3*8
	{
		TempMatrix[i] = new double[8];
	}
	
	EMatrix = new double *[3];
	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];
	
	BMatrix = new double *[3];		//3*8
	for (i=0; i<3; i++)
		BMatrix[i] = new double [8];
	
	StiffMat= new double*[8];
	for (i=0;i<8;i++)					// 8*8
		StiffMat[i] = new double[8];
	
	
	GaussB = new double **[4];		//3*8
	for (i=0;i<4;i++){
		GaussB[i] = new double *[3];
	}
	
	for (i=0;i<4;i++){
		for (j=0;j<3;j++){
			GaussB[i][j] = new double [8];
		}
	}
	
	
	
	StrainMatrix = new double *[4];
	for (i=0;i<4;i++) StrainMatrix[i] = new double[3];
	
	StressMatrix = new double *[4];
	for (i=0;i<4;i++) StressMatrix[i] = new double[3];
	
	
	
}
///////////////////////////////////////

void Plane4::CalculateMatrix(Node *DNodes, Material **DMat, long Loop)
{
	
	/////*****************************************////////
	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 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
	
	
	//////////////////*****************************/////////////////////	
	//Calculation of EMatrix
	EMatrixConst=DMat[Mat]->Ex[Loop]/(1-pow(DMat[Mat]->Nuxy[Loop],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[Loop]*EMatrixConst);
	EMatrix[1][0]=(DMat[Mat]->Nuxy[Loop]*EMatrixConst);
	EMatrix[1][1]=EMatrixConst;
	EMatrix[2][2]=(0.5*(1-DMat[Mat]->Nuxy[Loop])*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;
		
		
		
		
		for (j=0; j<3;j++) {
			for (k=0;k<8;k++){
				GaussB[i][j][k] = BMatrix[j][k];
			}
		}
		
	
		// Calculation of the Transpose of BMatrix[][]
		matrix_transpose(BtMatrix,BMatrix,8,3);
		
		matrix_multiply(EMatrix, BMatrix, TempMatrix,3,3,8);		//first and last values give rows and 
		matrix_multiply(BtMatrix, TempMatrix, StiffMat,8,3,8);		//columns respectively for resultant matrix
		
		// 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];
}



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

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");
}


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

void Plane4::CalculateStress(Node *DNodes, Gauss22 *DGauss22)
{
	int i,j,k;
	double *XStress;
	double *YStress;
	double *ShStress;
	double *XStrain;
	double *YStrain;
	double *ShStrain;

	
	XStress = new double[4];
	YStress = new double[4];
	ShStress = new double[4];
	XStrain = new double[4];
	YStrain = new double[4];
	ShStrain = new double[4];

	double *UMatrix = new double [8];
	
	for (j=0; j<4; j++){
		for (k=0; k<2; k++){
			UMatrix[2*j+k] = DNodes[NbLink[j]-1].u[k];
		}
	}
	
	
	//********* Shspe function to interpolate to node points  ( by converting factor sqrt3)********/
	
	
	//double *N;
	//N = new double [4];
	//double sq3;
	//sq3 = sqrt(3);
	//N[0] = 0.25*(1+sq3)*(1+sq3);
	//N[1] = 0.25*(1-sq3)*(1+sq3);
	//N[2] = 0.25*(1-sq3)*(1-sq3);
	//N[3] = 0.25*(1+sq3)*(1-sq3);
	
	double **Ng;
	Ng = new double *[4];
	for (i=0;i<4;i++){
		Ng[i] = new double [4];
	}
	
	
	for (i=0; i<4; i++){
		Ng[0][i] = DGauss22->N[i];
		Ng[3][i] = DGauss22->N[3-i];		//order of Ng4 is reverse of Ng1
	}
	

	Ng[1][0] = DGauss22->N[1];
	Ng[1][1] = DGauss22->N[0];
	Ng[1][2] = DGauss22->N[3];
	Ng[1][3] = DGauss22->N[2];
	
	Ng[2][0] = DGauss22->N[2];
	Ng[2][1] = DGauss22->N[3];
	Ng[2][2] = DGauss22->N[0];
	Ng[2][3] = DGauss22->N[1];
	

	//********  Calculate Strains at gauss points  ( epel = [B]*{d} ) *********/
	for (i=0;i<4;i++){
		pmatrix_vector_multiply(GaussB,UMatrix,3,8,StrainMatrix,4);
	}
	
	for (i=0;i<4;i++){
		XStrain[i] = StrainMatrix[i][0];	
		YStrain[i] = StrainMatrix[i][1];	
		ShStrain[i] = StrainMatrix[i][2];
	}
	
	for(i=0;i<4;i++){
		StrainMatrix[i][0] = Ng[i][0]*XStrain[0]+Ng[i][1]*XStrain[1]+Ng[i][2]*XStrain[2]+Ng[i][3]*XStrain[3];
		StrainMatrix[i][1] = Ng[i][0]*YStrain[0]+Ng[i][1]*YStrain[1]+Ng[i][2]*YStrain[2]+Ng[i][3]*YStrain[3];
		StrainMatrix[i][2] = Ng[i][0]*ShStrain[0]+Ng[i][1]*ShStrain[1]+Ng[i][2]*ShStrain[2]+Ng[i][3]*ShStrain[3];
	}

	
	matrix_pvector_multiply(EMatrix, StrainMatrix, 3,3, StressMatrix,4);
	
	for (i=0;i<4;i++){
		XStress[i] = StressMatrix[i][0];	
		YStress[i] = StressMatrix[i][1];	
		ShStress[i] = StressMatrix[i][2];
	}
	
	for(i=0;i<4;i++){
		StressMatrix[i][0] = Ng[i][0]*XStress[0]+Ng[i][1]*XStress[1]+Ng[i][2]*XStress[2]+Ng[i][3]*XStress[3];
		StressMatrix[i][1] = Ng[i][0]*YStress[0]+Ng[i][1]*YStress[1]+Ng[i][2]*YStress[2]+Ng[i][3]*YStress[3];
		StressMatrix[i][2] = Ng[i][0]*ShStress[0]+Ng[i][1]*ShStress[1]+Ng[i][2]*ShStress[2]+Ng[i][3]*ShStress[3];
	}
	
	
	
	for (i=0;i<3;i++){
		for (j=0;j<4;j++){
			DNodes[(NbLink[j]-1)].Stress[i] += StressMatrix[i][i];
			DNodes[(NbLink[j]-1)].Strain[i] += StrainMatrix[i][i];

		}
	}

	//delete unwanted memory
	delete []XStress;	XStress=0;
	delete []YStress;	YStress=0;
	delete []ShStress;	ShStress=0;
	delete []XStrain;	XStrain=0;
	delete []YStrain;	YStrain=0;
	delete []ShStrain;	ShStrain=0;
	delete []UMatrix;	UMatrix=0;
	for (i=0;i<4;i++) {
		delete []Ng[i];
		Ng[i]=0;
	}
}