function y = SpaceFrameElementMass(Rou,A,J,node_axis1,node_axis2)
x1=node_axis1(1);y1=node_axis1(2);z1=node_axis1(3);
x2=node_axis2(1);y2=node_axis2(2);z2=node_axis2(3);
L = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
w1 = 1/3;
w2 = 13/35;
w3 = J/(3*A);
w4 = 11*L/210;
w5 = L*L/105;
w6 = 9/70;
w7 = 13*L/420;
w8 = L*L/140;
m1=zeros(12,12);
m1(1,1)=w1;
m1(2,2)=w2;
m1(3,3)=w2;
m1(4,4)=w3;
m1(5,3)=-w4;
m1(5,5)=w5;
m1(6,2)=w4;
m1(6,6)=w5;
m1(7,1)=w1/2;
m1(7,7)=w1;
m1(8,2)=w6;
m1(8,6)=w7;
m1(8,8)=w2;
m1(9,3)=w6;
m1(9,5)=-w7;
m1(9,9)=w2;
m1(10,4)=w3/2;
m1(10,10)=w3;
m1(11,3)=w7;
m1(11,5)=-w8;
m1(11,9)=w4;
m1(11,11)=w5;
m1(12,2)=-w7;
m1(12,6)=-w8;
m1(12,8)=-w4;
m1(12,12)=w5;
m2=m1';
mprime = m1+m2;
if x1 == x2 & y1 == y2
   if z2 > z1
      Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
   else
      Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
   end
else
   CXx = (x2-x1)/L;
	CYx = (y2-y1)/L;
	CZx = (z2-z1)/L;
	D = sqrt(CXx*CXx + CYx*CYx);
	CXy = -CYx/D;
	CYy = CXx/D;
	CZy = 0;
	CXz = -CXx*CZx/D;
	CYz = -CYx*CZx/D;
	CZz = D;
	Lambda = [CXx CYx CZx ; CXy CYy CZy ; CXz CYz CZz];
end
R = [Lambda zeros(3) zeros(3) zeros(3) ; 
   zeros(3) Lambda zeros(3) zeros(3) ;
   zeros(3) zeros(3) Lambda zeros(3) ;
   zeros(3) zeros(3) zeros(3) Lambda];
y = R'*mprime*R;