//-------LINK1---------------------------------------------------------
// 2D-spar Element

class Link1					
{
private:
	long number;
	long NbLink[2];			// connectivity
	long Mat;				// Material
	long R;					// Real (for cross section area)
	double Angle;			// Angle between the x-axis and the bar
	double Length;			// lenght of the bar
	double Volume;
	double Mass;
	double StiffMat[4][4];
public:
	Link1(long num, long mat, long real, long* connectivity, Node *DNodes);
	void CalculateMatrix(Node *DNodes, Material **DMat, CrossSection **DCrossSection);
	void Link1::WriteStiffMatrix(FILE *file);
	void AssembleMatrix(double **Global,int MaxDOF,Node *DNodes);
};

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

Link1::Link1(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[1]-1)].ActiveDOF[0]=1;	// x DOF active node 2
	DNodes[(NbLink[1]-1)].ActiveDOF[1]=1;	// y DOF active
}

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

void Link1::CalculateMatrix(Node *DNodes, Material **DMat, CrossSection **DCrossSection)
{
	double C,S,K;
	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;
	K = DCrossSection[R]->Area*DMat[Mat]->Ex/Length;
	C=dx/Length;
	S=dy/Length;
	Angle = acos(C)*180/PI;
	StiffMat[0][0]= K*pow(C,2);
	StiffMat[0][1]= K*S*C;
	StiffMat[0][2]= -StiffMat[0][0];
	StiffMat[0][3]= -StiffMat[0][1];
	StiffMat[1][0]= StiffMat[0][1];
	StiffMat[1][1]= K*pow(S,2);
	StiffMat[1][2]= -StiffMat[0][1];
	StiffMat[1][3]= -StiffMat[1][1];
	StiffMat[2][0]= -StiffMat[0][0];
	StiffMat[2][1]= -StiffMat[0][1];
	StiffMat[2][2]= StiffMat[0][0];
	StiffMat[2][3]= StiffMat[0][1];
	StiffMat[3][0]= -StiffMat[0][1];
	StiffMat[3][1]= -StiffMat[1][1];
	StiffMat[3][2]= StiffMat[0][1];
	StiffMat[3][3]= StiffMat[1][1];
}

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

void Link1::AssembleMatrix(double **Global, int MaxDOF, Node *DNodes)
{
	int i,j,k,l;
	for(i=0;i<2;i++)			// step for column (matrix 2x2)
		for(j=0;j<2;j++)		// step for row (matrix 2x2)
			for(k=0;k<2;k++)	// step for line
				for(l=0;l<2;l++)
					Global[(DNodes[(NbLink[i]-1)].CumSum+k)][(DNodes[NbLink[j]-1].CumSum+l)]+=StiffMat[(i*2)+k][(j*2)+l];
}

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

void Link1::WriteStiffMatrix(FILE *file)
{
	int i,j;
	fprintf(file,"Element matrix %d, vol=%13.3E\n",number,Volume);
	for (i=0;i<4;i++)
		{
			for (j=0;j<4;j++)
			{
				fprintf(file,"%13.3E ",StiffMat[i][j]);
			}
			fprintf(file,"\n");
		}
		fprintf(file,"\n");
	}