// 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 **StiffMat;		// 
	double Area;		//Element's Area..
public:
	Plane3(long num, long mat, long real, long *connectivity, Node *DNodes);
	void CalculateMatrix(Node *DNodes, Material **DMat);
	void WriteStiffMatrix(FILE *file);
	void AssembleMatrix(double **Global,int MaxDOF,Node *DNodes);
};

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

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

	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
		}
}

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

void Plane3::CalculateMatrix(Node *DNodes, Material **DMat)
{
	int i,j;
	double DMatrixConst;
	double Det;

	double x1,x2,x3;
	double y1,y2,y3;

	double **DMatrix=new double*[3];
	double **BMatrix=new double*[3];		// B Matrix
	double **BtMatrix=new double*[6];		// Transpose of B Matrix
	double **TempMatrix=new double*[3];

	for (i=0;i<3;i++)			// 3*3
		DMatrix[i] = new double[3];

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

	StiffMat= new double*[6];
	
	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 D
	DMatrixConst=DMat[Mat]->Ex/(1-pow(DMat[Mat]->Nuxy,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*DMatrixConst;
	DMatrix[1][0]=DMat[Mat]->Nuxy*DMatrixConst;
	DMatrix[1][1]=DMatrixConst;
	DMatrix[2][2]=0.5*(1-DMat[Mat]->Nuxy)*DMatrixConst;

	//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];


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