//-------LINK1---------------------------------------------------------
// 2D-spar Element
#define PI 3.14159265359


class Link1					
{
private:
	long number;
	long NbLink[2];			// connectivity
	long Mat;				// Material
	long R;				// Real (for cross section area)
	double C,S;
	double dx,dy;


	double Angle;			// Angle between the x-axis and the bar
	double Length;			// length of the element
	double Volume;
	double Mass;
	double StiffMat[4][4];
	double Strain;
	double Stress;
public:
	Link1(long num, long mat, long real, long* connectivity, Node *DNodes, CrossSection **DCrossSection);
	void CalculateMatrix(Node *DNodes, Material **DMat, CrossSection **DCrossSection, long Loop);
	void Link1::WriteStiffMatrix(FILE *file);
	void AssembleMatrix(double **Global,int MaxDOF,Node *DNodes);
	void CalculateStrains(Node *DNodes, Material **DMat, long Loop);
	double GetStrain() const {return Strain; }
	double GetStress() const { return Stress; }
};

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

Link1::Link1(long num, long mat, long real, long* connectivity, Node *DNodes, CrossSection **DCrossSection)
{
	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
	
	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.");
	C=dx/Length;
	S=dy/Length;
	Angle = acos(C)*180/PI;
	Strain  =0.0;
	Stress = 0.0;
}

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

void Link1::CalculateMatrix(Node *DNodes, Material **DMat, CrossSection **DCrossSection, long Loop)
{
	double K;

	Mass=Volume*DMat[Mat]->Density[Loop];
	K = DCrossSection[R]->Area*DMat[Mat]->Ex[Loop]/Length;

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

/**************************************************/
void Link1::CalculateStrains(Node *DNodes, Material **DMat, long Loop)
{
	double duy, dux;
	double dLength;

	dux = DNodes[NbLink[1]-1].u[0]-DNodes[NbLink[0]-1].u[0];
	duy = DNodes[NbLink[1]-1].u[1]-DNodes[NbLink[0]-1].u[1];
	dLength = dux*C+duy*S;
	
	Strain = dLength/Length;
	Stress = Strain*DMat[Mat]->Ex[Loop];
}