// beam2.h
//
//--------Beam2-------------------
//

class Beam2
{
private:
	long number;
	long NbLink[2];			// connectivity
	long Mat;				// Material
	long R;					// Real
	double EL;				// EL
	double Volume;
	double Mass;
	double Angle;			// Angle between the x-axis and the beam
	double Length;
	double StiffMat[6][6];	// Stiffness Matrix dim = 3 (3*2Nodes)
public:
	Beam2(long num, long mat, long real, long *connectivity, Node *DNodes);
	void CalculateMatrix(Node *DNodes, Material **DMat, CrossSection **DCrossSection);
	void AssembleMatrix(double **Global,int MaxDOF,Node *DNodes);
	void WriteStiffMatrix(FILE *file);
};

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

Beam2::Beam2(long num,long mat, long real, long *connectivity, Node *DNodes)
{
	number = num;
	Mat = mat;
	R = real;

	NbLink[0] = connectivity[0];
	NbLink[1] = connectivity[1];

	DNodes[(NbLink[0]-1)].ActiveDOF[0]=1;	// x DOF active node 1
	DNodes[(NbLink[0]-1)].ActiveDOF[1]=1;	// y DOF active
	DNodes[(NbLink[0]-1)].ActiveDOF[5]=1;	// rotz DOF active
	DNodes[(NbLink[1]-1)].ActiveDOF[0]=1;	// x DOF active node 2
	DNodes[(NbLink[1]-1)].ActiveDOF[1]=1;	// y DOF active
	DNodes[(NbLink[1]-1)].ActiveDOF[5]=1;	// rotz DOF active
}

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

void Beam2::CalculateMatrix(Node *DNodes, Material **DMat, CrossSection **DCrossSection)
{
	int i,j;
	double A, B, K, M, F, G, P, H, Q,Z; // temporary variables
	double C,S;
	double dx,dy;
	dx=DNodes[(NbLink[1] - 1)].nx[0]-DNodes[(NbLink[0] - 1)].nx[0];
	dy=DNodes[(NbLink[1] - 1)].nx[1]-DNodes[(NbLink[0] - 1)].nx[1];
	Length=sqrt(pow(dx,2)+pow(dy,2));
	Volume=Length*DCrossSection[R]->Area;
	if(Volume<=0.0)
		nerror("Element volume must be greater than 0.");
	Mass=Volume*DMat[Mat]->Density;
	C=dx/Length;
	S=dy/Length;
	Angle = acos(C)*180/PI;
	EL=DMat[Mat]->Ex/Length;	// E/L ratio Young modulus over the length
	//Previous Calculations : Check them 'cos I'm not sure of the result...
	A=4*EL*DCrossSection[R]->Izz;				//A=Area Pense que ce n'est pas bon du tout;
	Z=DCrossSection[R]->Area*EL;	// Verify this is the area use instead of A
	B=2*EL*DCrossSection[R]->Izz; 
	K=12*EL*DCrossSection[R]->Izz/Length/Length;
	M = 6*EL*DCrossSection[R]->Izz/Length; 
	F= Z*pow(C,2)+ K*pow(S,2);
	G= (Z-K)*C*S;
	P = Z*pow(S,2)+ K*pow(C,2);
	H = -M*S;
	Q= M*C;
	StiffMat[0][0]= F;
	StiffMat[0][1]= G;
	StiffMat[0][2]= H;
	StiffMat[0][3]= -F;
	StiffMat[0][4]= -G;
	StiffMat[0][5]= H;
	StiffMat[1][1]= P;
	StiffMat[1][2]= Q;
	StiffMat[1][3]= -G;
	StiffMat[1][4]= -P;
	StiffMat[1][5]= Q;
	StiffMat[2][2]= A;
	StiffMat[2][3]= -H;
	StiffMat[2][4]= -Q;
	StiffMat[2][5]= B;		
	StiffMat[3][3]= F;
	StiffMat[3][4]= G;
	StiffMat[3][5]= -H;
	StiffMat[4][4]= P;
	StiffMat[4][5]= -Q;
	StiffMat[5][5]= A;
	// Copy matrix symmetric half
	for (i=0;i<6;i++)
		for (j=i+1;j<6;j++)
			StiffMat[j][i]=StiffMat[i][j];
}

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

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

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

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

void Beam2::WriteStiffMatrix(FILE *file)
{
	int i,j;
	fprintf(file,"Element matrix %d, vol=%13.3E\n",number,Volume);
	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");
}