% Model geometry
Coord (1,:) = [0  0  0];
n_fix = 1;

% Eléments de la poutre flexibles
n = size(Coord,1);  
ne = 0;
dL = L/np;
for i = 1:np
    Coord (n+i,:) = [L/np*i  0  0];
    Elem(ne+i,:) = [n+i-1 n+i];
end
n_extreme = np;

% Masselettes ...
for i = 1:np
    if (abs(L1 - Coord(n+i,1)) <= dL/2)
        Coord(n+i,1) = L1;
        Coord(n+i-1,1) = L1-LM/2;
        Coord(n+i+1,1) = L1+LM/2;
        n_1 = n+i;
        eMass = [n+i n+i+1];
        break
    end
end
    
for i = 1:np
    if (abs(L2 - Coord(n+i,1)) <= dL/2)
        Coord(n+i,1) = L2;
        Coord(n+i-1,1) = L2-LM/2;
        Coord(n+i+1,1) = L2+LM/2;
        n_M = n+i;
        eMass = [eMass  n+i n+i+1];
        break
    end
end

for i = 1:np
    if (abs(L3 - Coord(n+i,1)) <= dL/2)
        Coord(n+i,1) = L3;
        Coord(n+i-1,1) = L3-LM/2;
        Coord(n+i+1,1) = L3+LM/2;
        n_2 = n+i;
        eMass = [eMass  n+i n+i+1];
        break
    end
end

eMass
eFlex = setdiff([1:np], eMass)

% Eléments rigides pour visualisation
n = size(Coord,1);  ne = size(Elem,1);
Coord(n+1,:) = Coord(n_1,:) + [0 0 HM/2];
Coord(n+2,:) = Coord(n_1,:) - [0 0 HM/2];
Elem(ne+1,:) = [ n_1 n+1 ];
Elem(ne+2,:) = [ n_1 n+2 ];
eStiff = [ne+1  ne+2];

n = size(Coord,1);  ne = size(Elem,1);
Coord(n+1,:) = Coord(n_M,:) + [0 0 HM/2];
Coord(n+2,:) = Coord(n_M,:) - [0 0 HM/2];
Elem(ne+1,:) = [ n_M n+1 ];
Elem(ne+2,:) = [ n_M n+2 ];
eStiff = [eStiff  ne+1  ne+2];

n = size(Coord,1);  ne = size(Elem,1);
Coord(n+1,:) = Coord(n_2,:) + [0 0 HM/2];
Coord(n+2,:) = Coord(n_2,:) - [0 0 HM/2];
Elem(ne+1,:) = [ n_2 n+1 ];
Elem(ne+2,:) = [ n_2 n+2 ];
eStiff = [eStiff  ne+1  ne+2];

% Actionneur-1
n = size(Coord,1);  ne = size(Elem,1);
Coord(n+1,:) = Coord(n_1,:)+[0        -0.001        0];  
Coord(n+2,:) = Coord(n_1,:)+[0        -L_act+0.001  0];  
Coord(n+3,:) = Coord(n_1,:)+[0        -L_act        0];  
Coord(n+4,:) = Coord(n_1,:)+[L_motor  -L_act        0];
N_motor1 = n+4;
Elem (ne+1,:) = [n_1 n+1];
Elem (ne+2,:) = [n+1 n+2];
Elem (ne+3,:) = [n+2 n+3];
Elem (ne+4,:) = [n+3 n+4];
eAct  = [ne+2  ne+4];
eH    = [ne+1  ne+3];

% Actionneur-2
n = size(Coord,1);  ne = size(Elem,1);
Coord(n+1,:) = Coord(n_2,:)+[0        0.001        0];  
Coord(n+2,:) = Coord(n_2,:)+[0        L_act-0.001  0];  
Coord(n+3,:) = Coord(n_2,:)+[0        L_act        0];  
Coord(n+4,:) = Coord(n_2,:)+[L_motor  L_act        0];
N_motor2 = n+4;
Elem (ne+1,:) = [n_2 n+1];
Elem (ne+2,:) = [n+1 n+2];
Elem (ne+3,:) = [n+2 n+3];
Elem (ne+4,:) = [n+3 n+4];
eAct  = [eAct ne+2  ne+4];
eH    = [eH  ne+1  ne+3];