% FEM static run

% ------ Boundary conditions ---------------------------------------
bc = [];
[b,bc,nb] = fix_point (bc,n_fix,Dof);
[b,bc,nb] = fix_xyz   (bc,N_motor1,Dof);
[b,bc,nbc] = fix_1d (bc,N_motor1,Dof,[4 5]);
[b,bc,nb] = fix_xyz   (bc,N_motor2,Dof);
[b,bc,nbc] = fix_1d (bc,N_motor2,Dof,[4 5]);

%  ------ Load - Z moment ----------------------------------
p = zeros(size(K,1),1);
i = N_motor1
   dof_exc = (i-1)*n_dof+6;
   p(dof_exc) = 0.1;   % Nm

[X,R,xyzF] = fe_stat (K,p,b,n_dof,n_nodes);

disp ('Moment 0.1 Nm about Z');
disp ('xyz n_1 (mm)');
disp (xyzF(n_1,1:3)*1000)

disp ('xyz beam extreme (mm)');
disp (xyzF(n_extreme,1:3)*1000)

Edb = extract (Edof,X);	

figure; femdraw2 ([Coord(:,1) Coord(:,2)],Ex,Ey,[2 4 2 4 0.]); 
Edbxy = [Edb(:,1) Edb(:,2) Edb(:,6) Edb(:,7) Edb(:,8) Edb(:,12)];
femdisp2 (Ex,Ey,Edbxy,[1 5 0 5],1); ylabel('y');
title('Static analysis 0.1 Nm torque');