// Plane3.h

//-Plane3 ----------------------------------------------------------
// Constant Strain Triangle

class Plane3
{
private:
	long number;			
	long NbLink[3];			// 3 connections
	long Mat;				// Material
	long R;					// Real
	double Thickness;		//Thickness of the element
	double Area;			//Element's Area..
	double DMatrixConst;
	double Det;

	double x1,x2,x3;
	double y1,y2,y3;
	double **StiffMat;
	double **DMatrix;
	double **BMatrix;		// B Matrix
	double **BtMatrix;		// Transpose of B Matrix
	double **TempMatrix;

public:
	Plane3(long num, long mat, long real, long *connectivity, Node *DNodes);
	void CalculateMatrix(Node *DNodes, Material **DMat, long Loop);
	void WriteStiffMatrix(FILE *file);
	void AssembleMatrix(double **Global,int MaxDOF,Node *DNodes);
	void CalculateStrains(Node *DNodes);
	double GetStrain(int i) const { return Strain[i]; }
	double GetStress(int i) const { return Stress[i]; }
	double *Strain;
	double *Stress;
};

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

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

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

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

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

	DMatrix=new double*[3];	
	BMatrix=new double*[3];
	BtMatrix=new double*[6];
	TempMatrix=new double*[3];
	StiffMat= new double*[6];
	
	Strain = new double [3];
	for (i=0; i<3;i++){ 
		Strain[i] =0; }
	
	Stress = new double [3];
	for (i=0;i<3;i++){ 
		Stress[i] = 0; }
	
	for (i=0;i<3;i++)			// 3*3
		DMatrix[i] = new double[3];

	for (i=0;i<3;i++)			// 3*6
	{
		BMatrix[i] = new double[6];
		TempMatrix[i] = new double[6];
	}
	
	for (i=0;i<6;i++)			// 6*3
		BtMatrix[i] = new double[3];


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

	x1 = DNodes[NbLink[0]-1].nx[0];
	x2 = DNodes[NbLink[1]-1].nx[0];
	x3 = DNodes[NbLink[2]-1].nx[0];
	y1 = DNodes[NbLink[0]-1].nx[1];
	y2 = DNodes[NbLink[1]-1].nx[1];
	y3 = DNodes[NbLink[2]-1].nx[1];

	//Calculation of the determinant = 2 * Element Area:

	Det=(x2*y3-x3*y2)-(x1*y3-x3*y1)+(x1*y2-x2*y1);
	Area=Det/2.0;

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

	BMatrix[0][0]=(y2-y3)/Det; BMatrix[0][2]=(y3-y1)/Det; BMatrix[0][4]=(y1-y2)/Det;
	BMatrix[1][1]=(x3-x2)/Det; BMatrix[1][3]=(x1-x3)/Det; BMatrix[1][5]=(x2-x1)/Det;
	BMatrix[2][0]=BMatrix[1][1];
	BMatrix[2][2]=BMatrix[1][3];
	BMatrix[2][4]=BMatrix[1][5];
	BMatrix[2][1]=BMatrix[0][0];
	BMatrix[2][3]=BMatrix[0][2];
	BMatrix[2][5]=BMatrix[0][4];


}

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

void Plane3::CalculateMatrix(Node *DNodes, Material **DMat, long Loop)
{
	int i,j;

	// Calculation of D
	DMatrixConst=DMat[Mat]->Ex[Loop]/(1-pow(DMat[Mat]->Nuxy[Loop],2));

	for (i=0;i<3;i++)
		for (j=0;j<3;j++)
			DMatrix[i][j]=0;
	
	DMatrix[0][0]=DMatrixConst;
	DMatrix[0][1]=DMat[Mat]->Nuxy[Loop]*DMatrixConst;
	DMatrix[1][0]=DMat[Mat]->Nuxy[Loop]*DMatrixConst;
	DMatrix[1][1]=DMatrixConst;
	DMatrix[2][2]=0.5*(1-DMat[Mat]->Nuxy[Loop])*DMatrixConst;


	// Calculation of the Transpose of BMatrix[][]
	matrix_transpose(BtMatrix,BMatrix,6,3);
			
	matrix_multiply(DMatrix, BMatrix, TempMatrix,3,3,6);
	matrix_multiply(BtMatrix, TempMatrix, StiffMat,6,3,6);

	for (i=0;i<6;i++)
		for (j=0;j<6;j++)
			StiffMat[i][j]=StiffMat[i][j]*Thickness*Area;
}

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

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

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

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

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

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

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

void Plane3::CalculateStrains(Node *DNodes)
{
	int i;
	double *dMatrix;
	dMatrix = new double [6];
	dMatrix[0] = DNodes[NbLink[0]-1].u[0];
	dMatrix[1] = DNodes[NbLink[0]-1].u[1];
	dMatrix[2] = DNodes[NbLink[1]-1].u[0];
	dMatrix[3] = DNodes[NbLink[1]-1].u[1];
	dMatrix[4] = DNodes[NbLink[2]-1].u[0];
	dMatrix[5] = DNodes[NbLink[2]-1].u[1];

	matrix_vector_multiply (BMatrix, dMatrix, 3,6, Strain);
	matrix_vector_multiply (DMatrix, Strain, 3,3, Stress);
	delete []dMatrix;
	dMatrix =0;
	cout << "\n";
	for (i=0;i<3;i++){
		cout << Strain[i] << "\t" ;
	}
	cout << "\n";
	for (i=0;i<3;i++){
		cout << Stress[i] << "\t";
	}

}

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