%
%---------------------------------------------------------------------
%
%     MaLCoLM version 1.0, December 2019
%
%     Multi-spiral column simulation module
%     OOFEM extension module
%
%---------------------------------------------------------------------
%     Author:      Petr Havlasek
%     Institution: Department of Mechanics
%                  Faculty of Civil Engineering
%                  Czech Technical University in Prague
%---------------------------------------------------------------------
%     Acknowledgment:
% Financial support for this work was provided by the Technology
% Agency of the Czech Republic (TAČR), project number TF05000040
% (CeSTaR: Computer simulation and experimental validation -
% complex service for flexible and efficient design of pre-cast concrete
% columns with innovative multi-spiral reinforcement).
%
%
%

clear all
clf
cla

WRITE = 1;

%% DEFINITION OF GEOMETRY

prompt = {'Column width [m]','Column height [m]'};
dlg_title = 'Square column geometry';
num_lines = [1 length(dlg_title)+20];
default = {'0.6','0.075'};

options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';

dim = inputdlg(prompt,dlg_title,num_lines,default,options);

dim_X = str2double(dim{1});
dim_Y = dim_X;
dim_Z = str2double(dim{2});


%% ELEMENT SIZE DEFINITION - REGULAR AND STRUCTURE MESH ONLY

max_aspect_ratio = 2;

def_elem_size_horizontal = round(min([dim_X, dim_Y])/15 * 100.) / 100.;
def_elem_size_vertical = round( ( dim_Z)/4 * 100.) / 100.;

min_elem_size = min([def_elem_size_horizontal, def_elem_size_vertical]);

if (def_elem_size_horizontal > max_aspect_ratio * min_elem_size)
    def_elem_size_horizontal = max_aspect_ratio * min_elem_size;
end

if (def_elem_size_vertical > max_aspect_ratio * min_elem_size)
    def_elem_size_vertical = max_aspect_ratio * min_elem_size;
end


prompt = {'Horizontal direction [m]','Vertical direction [m]'};
dlg_title = 'Finite element size';
num_lines = [1 length(dlg_title)+20];
default = {num2str(def_elem_size_horizontal), num2str(def_elem_size_vertical)};

options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';

el_size = inputdlg(prompt,dlg_title,num_lines,default,options);

elem_size_X = str2double(el_size{1});
elem_size_Y = elem_size_X;
elem_size_Z = str2double(el_size{2});

min_elem_size = min([elem_size_X, elem_size_Y, elem_size_Z]);

if (elem_size_X > max_aspect_ratio * min_elem_size)
    elem_size_X = max_aspect_ratio * min_elem_size;
end
if (elem_size_Y > max_aspect_ratio * min_elem_size)
    elem_size_Y = max_aspect_ratio * min_elem_size;
end
if (elem_size_Z > max_aspect_ratio * min_elem_size)
    elem_size_Z = max_aspect_ratio * min_elem_size;
end

n_elem_X = ceil(dim_X/elem_size_X);
n_elem_Y = ceil(dim_Y/elem_size_Y);
n_elem_Z = ceil(dim_Z/elem_size_Z);


% new element size
elem_size_X = dim_X/n_elem_X;
elem_size_Y = dim_Y/n_elem_Y;
elem_size_Z = dim_Z/n_elem_Z;

x_coords = linspace(-dim_X/2,dim_X/2,n_elem_X+1);
y_coords = linspace(-dim_Y/2,dim_Y/2,n_elem_Y+1);
z_coords = linspace(0,dim_Z,n_elem_Z+1);

%% REINFORCEMENT TOPOLOGY

% estimate what the user wants to model
H = 75.e-3;
if (dim_Z <= 100.e-3)
    N_loops = 1.;
    H = dim_Z;
else
    N_loops = ceil(dim_Z/H);
    H = dim_Z/N_loops;
end

prompt = {'concrete cover [mm]','spiral pitch [mm]','diameter or small spiral [mm]'};
dlg_title = 'Reinforcement topology';
num_lines = [1 length(dlg_title)+20];
default = {'25.',num2str(H*1000.), '180.'};

options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';

el_size = inputdlg(prompt,dlg_title,num_lines,default,options);

cover = str2double(el_size{1})/1000.;
H = str2double(el_size{2})/1000.;
d_S = str2double(el_size{3})/1000.;

N_loops = ceil(dim_Z/H);
H = dim_Z/N_loops;

%% REINFORCEMENT BAR SIZES

[answer_S, answer_Ss, answer_L, answer_Ls, answer_V, answer_Vs] = rebar_dialog;

switch answer_S
    case '#3 (D10)'
        D_S = 10e-3;
        A_S = 71.E-6;
    case '#4 (D13)'
        D_S = 13e-3;
        A_S = 129.E-6;
    case '#5 (D16)'
        D_S = 16e-3;
        A_S = 200.E-6;
    case '#6 (D19)'
        D_S = 19e-3;
        A_S = 284.E-6;
    case '#7 (D22)'
        D_S = 22e-3;
        A_S = 387.E-6;
    case '#8 (D25)'
        D_S = 25e-3;
        A_S = 509.E-6;
    case '#9 (D29)'
        D_S = 29e-3;
        A_S = 645.E-6;
    case '#10 (D32)'
        D_S = 32e-3;
        A_S = 819.E-6;
end

switch answer_L
    case '#3 (D10)'
        D_L = 10e-3;
        A_L = 71.E-6;
    case '#4 (D13)'
        D_L = 13e-3;
        A_L = 129.E-6;
    case '#5 (D16)'
        D_L = 16e-3;
        A_L = 200.E-6;
    case '#6 (D19)'
        D_L = 19e-3;
        A_L = 284.E-6;
    case '#7 (D22)'
        D_L = 22e-3;
        A_L = 387.E-6;
    case '#8 (D25)'
        D_L = 25e-3;
        A_L = 509.E-6;
    case '#9 (D29)'
        D_L = 29e-3;
        A_L = 645.E-6;
    case '#10 (D32)'
        D_L = 32e-3;
        A_L = 819.E-6;
end

switch answer_V
    case '#3 (D10)'
        D_V = 10e-3;
        A_V = 71.E-6;
    case '#4 (D13)'
        D_V = 13e-3;
        A_V = 129.E-6;
    case '#5 (D16)'
        D_V = 16e-3;
        A_V = 200.E-6;
    case '#6 (D19)'
        D_V = 19e-3;
        A_V = 284.E-6;
    case '#7 (D22)'
        D_V = 22e-3;
        A_V = 387.E-6;
    case '#8 (D25)'
        D_V = 25e-3;
        A_V = 509.E-6;
    case '#9 (D29)'
        D_V = 29e-3;
        A_V = 645.E-6;
    case '#10 (D32)'
        D_V = 32e-3;
        A_V = 819.E-6;
end

switch answer_Ss
    case 'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%'
        steel_S_sig_0 = 473;
        steel_S_H = 7.5e3;
        steel_S_omega_c = 1.;
        steel_S_a = 6.;
        
    case 'SD420 fy 497 MPa, fu = 723 MPa, eps = 10%'
        steel_S_sig_0 = 479;
        steel_S_H = 8.2e3;
        steel_S_omega_c = 1.;
        steel_S_a = 6.;
        
    case 'B500A fy 500 MPa, fu = 525 MPa, eps = 2.5%'
        steel_S_sig_0 = 500;
        steel_S_H = 9.e3;
        steel_S_omega_c = 1.;
        steel_S_a = 13.;
    case 'B500B fy 500 MPa, fu = 540 MPa, eps = 5.0%'
        steel_S_sig_0 = 500;
        steel_S_H = 5.5e3;
        steel_S_omega_c = 1.;
        steel_S_a = 7.25;
    case 'B500C fy 500 MPa, fu = 575 MPa, eps = 7.5%'
        steel_S_sig_0 = 500;
        steel_S_H = 5.5e3;
        steel_S_omega_c = 1.;
        steel_S_a = 6.15;
    case 'fy 500 MPa, eps = infty'
        steel_S_sig_0 = 500;
        steel_S_H = 0.;
        steel_S_omega_c = 0.;
        steel_S_a = 0.;
end

switch answer_Ls
    case 'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%'
        steel_L_sig_0 = 473;
        steel_L_H = 7.5e3;
        steel_L_omega_c = 1.;
        steel_L_a = 6.;
        
    case 'SD420 fy 497 MPa, fu = 723 MPa, eps = 10%'
        steel_L_sig_0 = 479;
        steel_L_H = 8.2e3;
        steel_L_omega_c = 1.;
        steel_L_a = 6.;
        
    case 'B500A fy 500 MPa, fu = 525 MPa, eps = 2.5%'
        steel_L_sig_0 = 500;
        steel_L_H = 9.e3;
        steel_L_omega_c = 1.;
        steel_L_a = 13.;
    case 'B500B fy 500 MPa, fu = 540 MPa, eps = 5.0%'
        steel_L_sig_0 = 500;
        steel_L_H = 5.5e3;
        steel_L_omega_c = 1.;
        steel_L_a = 7.25;
    case 'B500C fy 500 MPa, fu = 575 MPa, eps = 7.5%'
        steel_L_sig_0 = 500;
        steel_L_H = 5.5e3;
        steel_L_omega_c = 1.;
        steel_L_a = 6.15;
    case 'fy 500 MPa, eps = infty'
        steel_L_sig_0 = 500;
        steel_L_H = 0.;
        steel_L_omega_c = 0.;
        steel_L_a = 0.;
end

switch answer_Vs
    case 'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%'
        steel_V_sig_0 = 473;
        steel_V_H = 7.5e3;
        steel_V_omega_c = 1.;
        steel_V_a = 6.;
        
    case 'SD420 fy 497 MPa, fu = 723 MPa, eps = 10%'
        steel_V_sig_0 = 479;
        steel_V_H = 8.2e3;
        steel_V_omega_c = 1.;
        steel_V_a = 6.;
        
    case 'B500A fy 500 MPa, fu = 525 MPa, eps = 2.5%'
        steel_V_sig_0 = 500;
        steel_V_H = 9.e3;
        steel_V_omega_c = 1.;
        steel_V_a = 13.;
    case 'B500B fy 500 MPa, fu = 540 MPa, eps = 5.0%'
        steel_V_sig_0 = 500;
        steel_V_H = 5.5e3;
        steel_V_omega_c = 1.;
        steel_V_a = 7.25;
    case 'B500C fy 500 MPa, fu = 575 MPa, eps = 7.5%'
        steel_V_sig_0 = 500;
        steel_V_H = 5.5e3;
        steel_V_omega_c = 1.;
        steel_V_a = 6.15;
    case 'fy 500 MPa, eps = infty'
        steel_V_sig_0 = 500;
        steel_V_H = 0.;
        steel_V_omega_c = 0.;
        steel_V_a = 0.;
end


%% CONCRETE

% mean compressive strength on cylinders [MPa]
fcm = 38;

prompt = {'f_{cm} [MPa]'};
dlg_title = 'Concrete compressive strength';
num_lines = [1 length(dlg_title)+30];
default = {num2str(fcm)};

options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';

strength = inputdlg(prompt,dlg_title,num_lines,default,options);

fcm = str2double(strength{1});

% characteristic value of compressive strength [MPa]
fck = fcm-8;
% mean tensile strength [MPa]
if (fck>50)
    fctm = 2.12*log(1+fcm/10);
else
    fctm = 0.3*fck^(2/3);
end
% fracture energy [N/m]
Gf = 73*fcm^0.18;
% initial modulus [GPa]
Eci = 21.5*(fcm/10)^(1/3);
% critical crack opening [m]
wf = 4.444*Gf/fctm*1e-6;
% Poisson's ratio
nu = 0.2;
% eccentricity
ecc = 0.525;
% initial value of hardening variable
kinit = 0.3;
% hardening modulus
Hp = 0.01;
% dilation factor
dilation = 0.85;
% hardening parameter Ahard
Ahard = 0.08;
% hardening parameter Ahard
Bhard = 0.003;
% hardening parameter Ahard
Chard = 2.;
% hardening parameter Ahard
Dhard = 1.e-6;
% softening parameter
Asoft = 15;
% softening parameter for compression
efc = 1.e-4;
% type of softening (1=bilinear)
% stype = 1;
% softening parameter for tension
wf1 = 0.15;
% softening parameter for tension
ft1 = 0.3;

prompt = {'Young''s modulus [GPa]', 'Poisson''s ratio [-]', 'f_t [MPa]', 'w_f [m]',...
    'eccentricity [-]', 'k_{init} initial value of hardening variable [-]', 'H_p hardening modulus [-]', 'dilation [-]',...
    'A_{hard} [-]', 'B_{hard} [-]', 'C_{hard} [-]', 'D_{hard} [-]',...
    'Asoft [-]', 'e_{fc} softening parameter for compression [-]', 'softening parameter for tension wf1 [-]', 'softening paramter for tension ft1 [-]' };
dlg_title = 'CDPM2 parameters';
num_lines = [1 length(dlg_title)+20];
default = { num2str(Eci), num2str(nu), num2str(fctm), num2str(wf),...
    num2str(ecc), num2str(kinit), num2str(Hp), num2str(dilation),...
    num2str(Ahard), num2str(Bhard), num2str(Chard), num2str(Dhard),...
    num2str(Asoft),num2str(efc),num2str(wf1),num2str(ft1)};

options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';

CDPM2 = inputdlg(prompt,dlg_title,num_lines,default,options);

% initial modulus [GPa]
Eci = str2double(CDPM2{1});
% Poisson's ratio
nu = str2double(CDPM2{2});
% tensile strength
fctm = str2double(CDPM2{3});
% critical crack opening [m]
wf = str2double(CDPM2{4});

% eccentricity
ecc = str2double(CDPM2{5});
% initial value of hardening variable
kinit = str2double(CDPM2{6});
% hardening modulus
Hp = str2double(CDPM2{7});
% dilation factor
dilation = str2double(CDPM2{8});

% hardening parameter Ahard
Ahard = str2double(CDPM2{9});
% hardening parameter Ahard
Bhard = str2double(CDPM2{10});
% hardening parameter Ahard
Chard = str2double(CDPM2{11});
% hardening parameter Ahard
Dhard = str2double(CDPM2{12});

% softening parameter
Asoft = str2double(CDPM2{13});
% softening parameter for compression
efc = str2double(CDPM2{14});
% softening parameter for tension
wf1 = str2double(CDPM2{15});
% softening parameter for tension
ft1 = str2double(CDPM2{16});


%% LOADING AND BOUNDARY CONDITIONS

[answer_BC, answer_LOAD] = loading_dialog;

% boundary conditions
BC_fully_fixed = 0;
BC_fixed = 0;
BC_periodic = 0;

LOAD_compression = 0;
LOAD_eccentric_compression = 0;
LOAD_bending = 0;
LOAD_compression_bending_bending = 0;

switch answer_BC
    case 'Fully fixed'
        BC_fully_fixed = 1;
    case 'Vertically fixed'
        BC_fixed = 1;
    case 'Periodic BC'
        BC_periodic = 1;
end

switch answer_LOAD
    case 'Axial compression (direct)'
        LOAD_compression = 1;
    case 'Eccentric compression (indirect)'
        LOAD_eccentric_compression = 1;
        
        if (BC_periodic)
            BC_periodic = 0;
            BC_fixed = 1;
        end
        
    case 'Pure bending (direct)'
        LOAD_bending = 1;
    case 'Compression and bending (direct)'
        LOAD_compression_bending_bending = 1;
end


%% PROPORTIONALITY OF LOADING

eccentric_X = 0.;
eccentric_Y = 0.;
if ( LOAD_eccentric_compression )
    
    prompt = {'x_F/B','y_F/B'};
    dlg_title = 'Vertical force position';
    num_lines = [1 length(dlg_title)+20];
    default = {num2str(1/6), '0.'};
    
    options.Resize='on';
    options.WindowStyle='normal';
    options.Interpreter='tex';
    
    loading = inputdlg(prompt,dlg_title,num_lines,default,options);
    
    eccentric_X = str2double(loading{1})*dim_X;
    eccentric_Y = str2double(loading{2})*dim_Y;
    
end

LOAD_eps = 1.;
LOAD_kappa1 = 0.;
LOAD_kappa2 = 0.;


if ( LOAD_bending)
    LOAD_eps = 1.;
end

if ( LOAD_bending)
    
    prompt = {'\kappa_1','\kappa_2'};
    dlg_title = 'Loading setup';
    num_lines = [1 length(dlg_title)+20];
    default = {'1.', '0.'};
    
    options.Resize='on';
    options.WindowStyle='normal';
    options.Interpreter='tex';
    
    loading = inputdlg(prompt,dlg_title,num_lines,default,options);
    
    LOAD_kappa1 = str2double(loading{1});
    LOAD_kappa2 = str2double(loading{2});
    
end


if ( LOAD_compression_bending_bending)
    
    prompt = {'\epsilon','\kappa_1','\kappa_2'};
    dlg_title = 'Loading setup';
    num_lines = [1 length(dlg_title)+20];
    default = {'1.','0.2', '0.'};
    
    options.Resize='on';
    options.WindowStyle='normal';
    options.Interpreter='tex';
    
    loading = inputdlg(prompt,dlg_title,num_lines,default,options);
    
    LOAD_eps = str2double(loading{1});
    LOAD_kappa1 = str2double(loading{2});
    LOAD_kappa2 = str2double(loading{3});
    
end



%%  PLOT ?

answer = questdlg('Plot geometry and FE mesh?', ...
    'Plot', ...
    'Yes','No','Reinforcement only','Reinforcement only');
% Handle response
switch answer
    case 'Yes'
        PLOT = 1.0;
    case 'No'
        PLOT = 0.0;
    case 'Reinforcement only'
        PLOT = 0.5;
end

PLOT_elements = floor(PLOT);
PLOT_bounding_box = ceil(PLOT);
PLOT_spiral = ceil(PLOT);


%% GENERATE NODES


% auxiliary variables - do not modify
count_all_nodes = 0;

n_nodes = (n_elem_X + 1) * (n_elem_Y + 1) * (n_elem_Z + 1);

nodes_count = 0;

% global NODES

NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).number = 0;
NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).coords = [0, 0, 0];

NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).BC_X = 0;
NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).BC_Y = 0;
NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).BC_Z = 0;

NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).master_X = 0;
NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).master_Y = 0;
NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).master_Z = 0;

NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).master_DOF_X = 0;
NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).master_DOF_Y = 0;
NODES((n_elem_X+1),(n_elem_Y+1),(n_elem_Z+1)).master_DOF_Z = 0;


for iz = 1: (n_elem_Z+1)
    for iy = 1: (n_elem_Y+1)
        for ix = 1: (n_elem_X+1)
            nodes_count = nodes_count + 1;
            NODES(ix,iy,iz).number = nodes_count;
            NODES(ix,iy,iz).coords = [(ix-1)*elem_size_X - dim_X/2, (iy-1)*elem_size_Y - dim_Y/2, (iz-1)*elem_size_Z];
        end
    end
end

count_all_nodes = nodes_count;

%% GENERATE BRICK ELEMENTS

n_elem_bricks = n_elem_X * n_elem_Y * n_elem_Z;

region_count = 1;
regionsets(region_count, :) = [1 n_elem_bricks];

count_all_elements = n_elem_bricks;


%% DRAW

if (PLOT_elements)
    
    figure(1)
    grey = 0.7 ;
    for i = 1:length(x_coords)
        plot3 (ones(1,5)*x_coords(i), ...
            [min(y_coords) max(y_coords) max(y_coords) min(y_coords) min(y_coords) ], ...
            [0, 0, dim_Z, dim_Z, 0], 'Color', [grey grey grey])
        hold on
    end
    
    for i = 1:length(y_coords)
        plot3 ( [min(x_coords) max(x_coords) max(x_coords) min(x_coords) min(x_coords) ], ...
            ones(1,5)*y_coords(i), ...
            [0, 0, dim_Z, dim_Z, 0], 'Color', [grey grey grey])
    end
    
    for i = 1:length(z_coords)
        plot3 ( [min(x_coords) max(x_coords) max(x_coords) min(x_coords) min(x_coords) ], ...
            [min(y_coords) min(y_coords) max(y_coords) max(y_coords) min(y_coords) ], ...
            ones(1,5)*z_coords(i), 'Color', [grey grey grey])
    end
    
    
end

if (PLOT_bounding_box)
    
    figure(1)
    
    plot3 (ones(1,5)*min(x_coords), ...
        [min(y_coords) max(y_coords) max(y_coords) min(y_coords) min(y_coords) ], ...
        [0, 0, dim_Z, dim_Z, 0], '-k')
    hold on
    axis equal
    
    plot3 (ones(1,5)*max(x_coords), ...
        [min(y_coords) max(y_coords) max(y_coords) min(y_coords) min(y_coords) ], ...
        [0, 0, dim_Z, dim_Z, 0], '-k')
    
    plot3 ( [min(x_coords) max(x_coords) max(x_coords) min(x_coords) min(x_coords) ], ...
        ones(1,5)*min(y_coords), ...
        [0, 0, dim_Z, dim_Z, 0], '-k')
    
    
    plot3 ( [min(x_coords) max(x_coords) max(x_coords) min(x_coords) min(x_coords) ], ...
        ones(1,5)*max(y_coords), ...
        [0, 0, dim_Z, dim_Z, 0], '-k')
    
end



%% definition of spiral reinforcement

% spiral radius (axial)
r_S = d_S/2 - D_S;
r_L = ( min(dim_X,dim_Y) - 2*cover - D_L)/2;

% center of small spirals wrt origin
DD_X = dim_X/2 - cover - r_S - D_S/2;
DD_Y = dim_Y/2 - cover - r_S - D_S/2;

nSpiral = 5;

nHanging = count_all_nodes;

V_transversal_reinf = 0.;

for spiral = 1:nSpiral
    
    if (spiral == 1)
        r = r_L;
        x0 = 0.0;
        y0 = 0.0;
        A_transversal = A_L;
    elseif (spiral == 2)
        r = r_S;
        x0 = -DD_X;
        y0 = -DD_Y;
        A_transversal = A_S;
    elseif (spiral == 3)
        r = r_S;
        x0 = DD_X;
        y0 = -DD_Y;
        A_transversal = A_S;
    elseif (spiral == 4)
        r = r_S;
        x0 = DD_X;
        y0 = DD_Y;
        A_transversal = A_S;
    elseif (spiral == 5)
        r = r_S;
        x0 = -DD_X;
        y0 = DD_Y;
        A_transversal = A_S;
    end
    
    % length of one pitch
    L_pitch = sqrt(H^2 + (2*pi*r)^2 );
    
    V_transversal_reinf = V_transversal_reinf + N_loops * L_pitch * A_transversal;
    
    % maximum length of 1 truss element
    %     L_max = min(dim_X, dim_Y)/20;
    % corresponding maximum angle
    %     theta_max = 2*pi*L_max/L_pitch;
    
    % better approach: reflects different diameters
    % 6*4 = 24 (one quarter divided into 6 segments)
    theta_max = 2*pi/24;
    
    % minimum length of 1 truss element
    L_min = 1.e-3;
    % corresponding minimum angle
    theta_min = 2*pi*L_min/L_pitch;
    
    
    x_intercept = [];
    for i = 1:length(x_coords)
        if ( (x_coords(i) >= x0-r) && (x_coords(i) <= x0+r) )
            x_intercept = [x_intercept x_coords(i)];
        end
    end
    
    y_intercept = [];
    for i = 1:length(y_coords)
        if ( (y_coords(i) >= y0-r) && (y_coords(i) <= y0+r) )
            y_intercept = [y_intercept y_coords(i)];
        end
    end
    
    z_intercept = [];
    for i = 1:length(z_coords)
        if ( (z_coords(i) >= 0) && (z_coords(i) <= H) )
            z_intercept = [z_intercept z_coords(i)];
        end
    end
    z_intercept = unique([z_intercept H]);
    
    theta = [];
    
    for i = 1:length(x_intercept)
        theta_1 = acos( (x_intercept(i)-x0)/r );
        theta_2 = pi + (pi-theta_1);
        % %         check:
        %          x_1 = x0 + r*cos(theta_1)
        %          x_2 = x0 + r*cos(theta_2)
        theta = [theta theta_1 theta_2];
    end
    
    for i = 1:length(y_intercept)
        theta_1 = asin( (y_intercept(i)-y0)/r );
        if ( theta_1 < 0 )
            theta_1 = theta_1 + 2*pi;
            theta_2 = pi + (2*pi - theta_1);
            
        elseif (theta_1 < pi/2)
            theta_2 = pi - theta_1;
            
        elseif (theta_1 < pi)
            theta_2 = pi-theta_1;
            
        elseif (theta_1 < 3*pi/2)
            theta_2 = 2*pi - (theta_1 - pi);
            
        else
            theta_2 = pi + (2*pi - theta_1);
            
        end
        
        %         check:
        %          y_1 = y0 + r*sin(theta_1)
        %          y_2 = y0 + r*sin(theta_2)
        theta = [theta theta_1 theta_2];
    end
    
    for i = 1:length(z_intercept)
        theta = [theta z_intercept(i)*2*pi/H];
    end
    
    % remove duplicates and sort
    theta = unique(theta);
    
    % check that all elements are shorter than L_max
    theta_supplement = [];
    for i = 1:(length(theta)-1)
        d_theta = theta(i+1)-theta(i);
        
        if (d_theta > theta_max)
            N = ceil(d_theta/theta_max); %number of segments current segment should be divided
            refined = linspace(theta(i), theta(i+1), N+1);
            refined(1) = [];
            refined(length(refined)) = [];
            theta_supplement = [theta_supplement  refined];
        end
    end
    
    if ~(isempty(theta_supplement))
        theta = [theta theta_supplement];
    end
    
    if (N_loops > 1)
        theta_orig = theta;
        for i = 1:(N_loops-1)
            theta = [theta theta_orig + i*2*pi*ones(1,length(theta_orig))];
        end
    end
    
    theta = unique(theta);
    
    % another check to guarantee that truss elements do not have zero
    % length
    i = 1;
    while( (i+1) < length(theta) )
        if ((theta(i+1) - theta(i)) < theta_min)
            theta(i+1)=[];
        else
            i = i+1;
        end
    end
    if ((theta(i+1) - theta(i)) < theta_min)
        theta(i) = [];
    end
    
    minimum = 1000.;
    for (i=1:length(theta)-1)
        if(theta(i+1)-theta(i))<minimum
            minimum = (theta(i+1)-theta(i));
            huhu = i;
        end
        
    end
    
    
    number = zeros(1,length(theta));
    X = zeros(1,length(theta));
    Y = zeros(1,length(theta));
    Z = zeros(1,length(theta));
    
    for i = 1:length(theta)
        nHanging = nHanging + 1;
        number(i) = nHanging;
        X(i) =  x0 + r * cos(theta(i));
        Y(i) =  y0 + r * sin(theta(i));
        Z(i) =  H * theta(i)/(2*pi);
    end
    SPIRAL(spiral).nr = number;
    SPIRAL(spiral).X = X;
    SPIRAL(spiral).Y = Y;
    SPIRAL(spiral).Z = Z;
    
    % update
    count_all_elements = count_all_elements + length(number) - 1;
    
    if (PLOT_spiral)
        figure(1)
        if (spiral == 1)
            plot3 (SPIRAL(spiral).X, SPIRAL(spiral).Y, SPIRAL(spiral).Z, '-b', 'LineWidth', 3)
        else
            plot3 (SPIRAL(spiral).X, SPIRAL(spiral).Y, SPIRAL(spiral).Z, '-r', 'LineWidth', 3)
        end
        
        hold on
    end
    
end

%%

% here r_L and r_S are already net radii (center to axis of rebar)

% auxiliary variable - distance of 2 vertical rebars from the center of the column
r2 = r_L + D_L/2 + D_V/2;

% auxiliary variable - distance of all vertical rebars from the center of
% the small spiral
r1 = r_S - D_S/2 - D_V/2;

% vertical spacing is the same for all rebars = length of the element
dL = dim_Z/n_elem_Z;

nRebar = 16;

%TODO - currently only square columns are supported
DD = DD_X;

% rebars on the diagonal
coord_1_x = DD + r1/sqrt(2);
coord_1_y = DD + r1/sqrt(2);

coord_2_x = DD - r1/sqrt(2);
coord_2_y = DD - r1/sqrt(2);


alpha = acos ( (2*DD^2 + r2^2 - r1^2 )/ (2*sqrt(2) * DD * r2) );
beta = pi/4-alpha;

coord_3_x = r2*sin(beta);
coord_3_y = r2*cos(beta);

coord_4_x = coord_3_y;
coord_4_y = coord_3_x;

% pre-allocate
aux = zeros(1,length(z_coords));
REBAR(nRebar).nr = aux;
REBAR(nRebar).X = aux;
REBAR(nRebar).Y = aux;
REBAR(nRebar).Z = aux;

for rebar = 1:nRebar
    if (rebar == 1)
        x0 = coord_1_x;
        y0 = coord_1_y;
    elseif (rebar == 2)
        x0 = coord_2_x;
        y0 = coord_2_y;
    elseif (rebar == 3)
        x0 = coord_3_x;
        y0 = coord_3_y;
    elseif (rebar == 4)
        x0 = coord_4_x;
        y0 = coord_4_y;
    elseif (rebar == 5)
        x0 = -coord_1_x;
        y0 = -coord_1_y;
    elseif (rebar == 6)
        x0 = -coord_2_x;
        y0 = -coord_2_y;
    elseif (rebar == 7)
        x0 = -coord_3_x;
        y0 = -coord_3_y;
    elseif (rebar == 8)
        x0 = -coord_4_x;
        y0 = -coord_4_y;
    elseif (rebar == 9)
        x0 = -coord_1_x;
        y0 = coord_1_y;
    elseif (rebar == 10)
        x0 = -coord_2_x;
        y0 = coord_2_y;
    elseif (rebar == 11)
        x0 = -coord_3_x;
        y0 = coord_3_y;
    elseif (rebar == 12)
        x0 = -coord_4_x;
        y0 = coord_4_y;
    elseif (rebar == 13)
        x0 = coord_1_x;
        y0 = -coord_1_y;
    elseif (rebar == 14)
        x0 = coord_2_x;
        y0 = -coord_2_y;
    elseif (rebar == 15)
        x0 = coord_3_x;
        y0 = -coord_3_y;
    elseif (rebar == 16)
        x0 = coord_4_x;
        y0 = -coord_4_y;
    end
    
    number = zeros(1,length(z_coords));
    for i = 1:length(z_coords)
        nHanging = nHanging + 1;
        number(i) = nHanging;
    end
    
    REBAR(rebar).nr = number;
    REBAR(rebar).X = x0 * ones(1,length(z_coords));
    REBAR(rebar).Y = y0 * ones(1,length(z_coords));
    REBAR(rebar).Z = z_coords;
    
    % update
    count_all_elements = count_all_elements + length(number) - 1;
    
    if (PLOT_spiral)
        figure(1)
        plot3 (REBAR(rebar).X, REBAR(rebar).Y, REBAR(rebar).Z, '-k', 'LineWidth', 4)
        hold on
    end
    
end


%% WRITE OUT

if (WRITE)
    % HEADER
    fid = fopen('multi-spiral_column.in', 'wt');
    fprintf(fid, 'multi-spiral_column.out\n');
    fprintf(fid, 'Simulation of multispiral rectangular column with dimensions (B1 x B2 x H) %3.3f x %3.3f x %3.3f m, cover = %3.3f m, pitch = %3.3f m, d_s = %3.3f m, D_L = %3.1f mm, D_S = %3.1f mm, D_V = %3.1f mm\n',...
        [dim_X dim_Y dim_Z cover H 2*r_S+D_S D_L*1000. D_S*1000. D_V*1000. ] );
    
    
    fprintf(fid, '# Reinforcement ratio VERTICAL DIRECTION RHO_V %3.3f %%\n',...
        [100. * A_V*nRebar / (dim_X * dim_Y) ] );
    
    fprintf(fid, '# Reinforcement ratio LATERAL DIRECTION RHO_L %3.3f %%\n',...
        [100. * V_transversal_reinf / (dim_X * dim_Y * dim_Z) ] );
    
    fprintf(fid, '#\n');
    fprintf(fid, '#\n');
    
    
    fprintf(fid, 'NonLinearStatic nmsteps 3 nsteps 1 smtype 8 nmodules 3 renumber 1\n');
    
    if ( ~LOAD_eccentric_compression )
        
        fprintf(fid, '#\n');
        fprintf(fid, 'nsteps 100 deltaT 1.0  rtolv 1.e-4 MaxIter 200 controlmode 1 stiffMode 2 lstype 4 donotfixload\n');
        fprintf(fid, 'nsteps 100 deltaT 1.0  rtolv 2.e-4 MaxIter 150 controlmode 1 stiffMode 2 lstype 4 donotfixload\n');
        fprintf(fid, 'nsteps 100 deltaT 1.0  rtolv 5.e-4 MaxIter 100 controlmode 1 stiffMode 2 lstype 4 donotfixload\n');
        fprintf(fid, '#\n');
        fprintf(fid, '#\n');
    else
        
        fprintf(fid, 'nsteps 100 controlmode 0 stepLength %3.3d minStepLength %3.3d rtolv 1e-4 hpc 2 1000001 3 hpcw 1 -1. Psi 0.0 MaxIter 200 stiffmode 2 reqiterations 20 lstype 4 initialguess 1 donotfixload\n', [(dim_Z*0.004)/100 (dim_Z*0.002)/100.]);
        fprintf(fid, 'nsteps 100 controlmode 0 stepLength %3.3d minStepLength %3.3d rtolv 2e-4 hpc 2 1000001 3 hpcw 1 -1. Psi 0.0 MaxIter 150 stiffmode 2 reqiterations 20 lstype 4 initialguess 1 donotfixload\n', [(dim_Z*0.004)/100 (dim_Z*0.002)/100.]);
        fprintf(fid, 'nsteps 100 controlmode 0 stepLength %3.3d minStepLength %3.3d rtolv 5e-4 hpc 2 1000001 3 hpcw 1 -1. Psi 0.0 MaxIter 100 stiffmode 2 reqiterations 20 lstype 4 initialguess 1 donotfixload\n', [(dim_Z*0.004)/100 (dim_Z*0.002)/100.]);
        
        fprintf(fid, '#\n');
        fprintf(fid, '#\n');
    end
    
    
    fprintf(fid, '# concrete\n');
    fprintf(fid, 'vtkxml tstep_step 5 domain_all vars 2 1 4 cellvars 2 27 13 primvars 1 1 stype 1 regionsets 1 1\n');
    fprintf(fid, '# spirals\n');
    fprintf(fid, 'vtkxml tstep_step 5 domain_all vars 2 1 4 cellvars 1 27 primvars 1 1 stype 1 regionsets 2 2 3\n');
    fprintf(fid, '# rebars\n');
    fprintf(fid, 'vtkxml tstep_step 5 domain_all vars 2 1 4 cellvars 1 27 primvars 1 1 stype 1 regionsets 1 4\n');
    fprintf(fid, '#\n');
    fprintf(fid, '#\n');
    
    fprintf(fid, 'domain 3D\n');
    fprintf(fid, 'OutputManager tstep_all dofman_output {1000000 1000001 1000002 1000003}\n');
    fprintf(fid, '#\n');
    fprintf(fid, 'ndofman %3d nelem %3d ncrosssect 4 nmat 4 nbc 5 nic 0 nltf 5 nset 4\n', [nHanging+4 count_all_elements] );
    fprintf(fid, '#\n');
    
    % NODES
    % nodes - basic mesh
    
    for iz = 1: (n_elem_Z+1)
        for iy = 1: (n_elem_Y+1)
            for ix = 1: (n_elem_X+1)
                
                if (iz ~= (n_elem_Z+1) )
                    % standard node
                    fprintf(fid, 'node %3d coords 3 %2.6f %2.6f %2.6f', ...
                        [NODES(ix,iy,iz).number NODES(ix,iy,iz).coords] );
                    
                    % boundary conditions if applicable
                    if (iz == 1)
                        if (BC_fully_fixed)
                            fprintf(fid, ' bc 3 1 1 1');
                            
                        elseif (BC_fixed)
                            if ( (ix == 1) && (iy == 1) )
                                fprintf(fid, ' bc 3 1 1 1');
                            elseif ( (ix == (n_elem_X+1) ) && (iy == 1) )
                                fprintf(fid, ' bc 3 0 1 1');
                            else
                                fprintf(fid, ' bc 3 0 0 1');
                            end
                            
                        elseif (BC_periodic)
                            if ( (ix == 1) && (iy == 1) )
                                fprintf(fid, ' bc 3 1 1 1');
                            elseif ( (ix == (n_elem_X+1) ) && (iy == 1) )
                                fprintf(fid, ' bc 3 0 1 0');
                            end
                        end
                    end
                    
                    % end of line
                    fprintf(fid, '\n')
                    
                    
                else
                    % nodes on top surface - link to bottom surface and
                    % supplementary 3 nodes
                    fprintf(fid, 'slavenode %3d coords 3 %2.6f %2.6f %2.6f doftype 3 2 2 2 masterDofMan 4 %3d 1000001 1000002 1000003 weights 4 1. 1. %2.6f %2.6f\n', ...
                        [NODES(ix,iy,iz).number NODES(ix,iy,iz).coords NODES(ix,iy,1).number NODES(ix,iy,iz).coords(1)/(dim_X/2) NODES(ix,iy,iz).coords(2)/(dim_Y/2) ] );
                end
                
            end
        end
    end
    
    
    % nodes -  hanging nodes - spirals
    fprintf(fid, '#\n#\n');
    fprintf(fid, '# HANGING NODES - SPIRALS\n');
    fprintf(fid, '#\n#\n');
    
    for spiral = 1:nSpiral
        if spiral == 1
            fprintf(fid, '#\n# large spiral (%1d)\n#\n', spiral);
        else
            fprintf(fid, '#\n# small spiral (%1d)\n#\n', spiral);
        end
        for i = 1:length(SPIRAL(spiral).nr)
            fprintf(fid, 'hangingNode %3d coords 3 %2.8f %2.8f %2.8f dofType 3 2 2 2\n', ...
                [SPIRAL(spiral).nr(i) SPIRAL(spiral).X(i) SPIRAL(spiral).Y(i) SPIRAL(spiral).Z(i)]);
        end
    end
    
    % nodes -  hanging nodes - longitudinal
    fprintf(fid, '#\n#\n');
    fprintf(fid, '# HANGING NODES - VERTICAL BARS\n');
    fprintf(fid, '#\n#\n');
    
    for rebar = 1:nRebar
        for i = 1:length(REBAR(rebar).nr)
            fprintf(fid, 'hangingNode %3d coords 3 %2.6f %2.6f %2.6f dofType 3 2 2 2\n', ...
                [REBAR(rebar).nr(i) REBAR(rebar).X(i) REBAR(rebar).Y(i) REBAR(rebar).Z(i)]);
        end
    end
    
    % nodes -  supplementary nodes for loading
    
    % load for indirect control
    fprintf(fid, '#\n#\n');
    fprintf(fid, 'hangingNode %3d coords 3 %2.6f %2.6f %2.6f dofType 3 2 2 2', ...
        [1000000 eccentric_X eccentric_Y dim_Z]);
    if (LOAD_eccentric_compression)
        fprintf(fid, ' load 1 5');
    end
    fprintf(fid, '\n');
    fprintf(fid, '#\n#\n');
    
    % dummy node - axial compression - loading controlled by ltf 2
    fprintf(fid, '# dummy node - compression\n');
    fprintf(fid, 'node %3d coords 3 %2.6f %2.6f %2.6f', [1000001 0. 0. dim_Z]);
    if (LOAD_compression)
        fprintf(fid, ' bc 3 1 1 2');
        
    elseif (LOAD_eccentric_compression)
        fprintf(fid, ' bc 3 0 0 0');
        
    elseif (LOAD_bending)
        fprintf(fid, ' bc 3 0 0 0');
        
    elseif (LOAD_compression_bending_bending)
        fprintf(fid, ' bc 3 1 1 2');
    end
    fprintf(fid, '\n');
    
    % dummy node bending about y-axis - loading controlled by ltf 3
    fprintf(fid, '# dummy node - bending about y-axis\n');
    fprintf(fid, 'node %3d coords 3 %2.6f %2.6f %2.6f', [1000002 dim_X/2 0. dim_Z]);
    if (LOAD_compression) % i.e. no rotation
        fprintf(fid, ' bc 3 1 1 1');
        
    elseif (LOAD_eccentric_compression)
        fprintf(fid, ' bc 3 0 0 0');
        
    elseif (LOAD_bending)
        fprintf(fid, ' bc 3 1 1 3');
        
    elseif (LOAD_compression_bending_bending)
        fprintf(fid, ' bc 3 1 1 3');
    end
    fprintf(fid, '\n');
    
    % dummy node bending about x-axis - loading controlled by ltf 4
    fprintf(fid, '# dummy node - bending about x-axis\n');
    fprintf(fid, 'node %3d coords 3 %2.6f %2.6f %2.6f', [1000003 0. dim_Y/2 dim_Z]);
    if (LOAD_compression) % i.e. no rotation
        fprintf(fid, ' bc 3 1 1 1');
        
    elseif (LOAD_eccentric_compression)
        fprintf(fid, ' bc 3 0 0 0');
        
    elseif (LOAD_bending)
        fprintf(fid, ' bc 3 1 1 4');
        
    elseif (LOAD_compression_bending_bending)
        fprintf(fid, ' bc 3 1 1 4');
    end
    fprintf(fid, '\n');
    
end



% ELEMENTS

% elements - basic structured mesh

fprintf(fid, '#\n#\n#\n# Concrete\n#\n#\n');
elem_count = 0;

for iz = 1: (n_elem_Z)
    for iy = 1: (n_elem_Y)
        for ix = 1: (n_elem_X)
            
            elem_count = elem_count + 1;
            
            ELEM_NODES = [NODES(ix, iy ,iz+1).number, NODES(ix, iy+1 ,iz+1).number, NODES(ix+1, iy+1 ,iz+1).number, NODES(ix+1, iy ,iz+1).number,...
                NODES(ix, iy ,iz).number, NODES(ix, iy+1 ,iz).number, NODES(ix+1, iy+1 ,iz).number, NODES(ix+1, iy ,iz).number];
            
            fprintf(fid, 'LSpaceBB %3d   nodes 8 %4d %4d %4d %4d %4d %4d %4d %4d\n', ...
                [elem_count ELEM_NODES]);
            
        end
    end
end

% elements - spirals
aux_elem_number = elem_count;

fprintf(fid, '#\n#\n#\n# Spirals\n#\n#\n');
for spiral = 1:nSpiral
    if spiral == 1
        fprintf(fid, '#\n# large spiral\n#\n');
    else
        fprintf(fid, '#\n# small spiral (%1d)\n#\n', spiral);
    end
    for i = 1:(length(SPIRAL(spiral).nr) - 1)
        aux_elem_number = aux_elem_number+1;
        fprintf(fid, 'Truss3d %3d nodes 2 %4d %4d\n', ...
            [aux_elem_number SPIRAL(spiral).nr(i) SPIRAL(spiral).nr(i+1)]);
    end
end

% elements - vertical rebars
fprintf(fid, '#\n#\n#\n# Rebars\n#\n#\n');
for rebar = 1:nRebar
    for i = 1:(length(REBAR(rebar).nr) - 1)
        aux_elem_number = aux_elem_number+1;
        fprintf(fid, 'Truss3d %3d nodes 2 %4d %4d\n', ...
            [aux_elem_number REBAR(rebar).nr(i) REBAR(rebar).nr(i+1)]);
    end
end

% CROSS SECTIONS

fprintf(fid, '#\n#\n#\n# CROSS-SECTIONS\n#\n#\n');
fprintf(fid, '#\n# concrete\n#\n');
fprintf(fid, 'SimpleCS 1 material 1 set 1\n');
fprintf(fid, '#\n# large spiral\n#\n');
fprintf(fid, 'SimpleCS 2 area %6.6d material 2 set 2\n', A_L);
fprintf(fid, '#\n# small spiral\n#\n');
fprintf(fid, 'SimpleCS 3 area %6.6d material 3 set 3\n', A_S);
fprintf(fid, '#\n# vertical rebars\n#\n');
fprintf(fid, 'SimpleCS 4 area %6.6d material 4 set 4\n', A_V);

% MATERIALS

fprintf(fid, '#\n#\n#\n# MATERIALS\n#\n#\n');

fprintf(fid, '# CONCRETE_DEFINITION\n');
%fprintf(fid, 'Con2DPM 1 d 0 E 30.30e3 n 0.2 talpha 0. fc 28. ft 2.21 hp 0.01 asoft 15. wf 268.e-6\n');

fprintf(fid, 'Con2DPM 1 d 0 E %3.3f n %3.3f talpha 0. fc %3.3f ft %3.3f wf %3.3d ecc %1.3f kinit %1.3f Hp %3.3f dilation %1.3f Ahard %3.3d Bhard %3.3d Chard %3.3f Dhard %3.3f Asoft %3.3f efc %3.3d stype 1 wf1 %3.3f ft1 %3.3f\n',...
    [Eci*1000. nu fcm fctm wf ecc kinit Hp dilation Ahard Bhard Chard Dhard Asoft efc wf1 ft1]);
fprintf(fid, '#\n');
fprintf(fid, '# large SPIRAL_DEFINITION - parameters from Kuo''s master thesis (average for transversal reinforcement)\n');
fprintf(fid, 'MisesMat 2 d 0. talpha 0. E 200.e3 n 0.3 sig0 %3.3d H %3.3f omega_crit %3.3f a %3.3f\n',...
    [steel_L_sig_0 steel_L_H steel_L_omega_c steel_L_a] );
fprintf(fid, '# small SPIRAL_DEFINITION - parameters from Kuo''s master thesis (average for transversal reinforcement)\n');
fprintf(fid, 'MisesMat 3 d 0. talpha 0. E 200.e3 n 0.3 sig0 %3.3d H %3.3f omega_crit %3.3f a %3.3f\n',...
    [steel_S_sig_0 steel_S_H steel_S_omega_c steel_S_a] );
fprintf(fid, '#\n');
fprintf(fid, '# REBAR_DEFINITION\n');
% fprintf(fid, 'MisesMat 4 d 0. talpha 0. E 200.e3 n 0.3 sig0 497. H 8.2e3 omega_crit 1. a 6.\n');
fprintf(fid, 'MisesMat 4 d 0. talpha 0. E 200.e3 n 0.3 sig0 %3.3d H %3.3f omega_crit %3.3f a %3.3f\n',...
    [steel_V_sig_0 steel_V_H steel_V_omega_c steel_V_a] );

% BOUNDARY CONDITIONS
fprintf(fid, '#\n#\n#\n# BOUNDARY CONDITIONS\n#\n#\n');
fprintf(fid, 'BoundaryCondition  1 loadTimeFunction 1 prescribedvalue 0.\n');
fprintf(fid, '#\n');
% eps = 0.01 = w / dim_Z
fprintf(fid, 'BoundaryCondition  2 loadTimeFunction 2 prescribedvalue %3.4d\n', -0.01*dim_Z);
fprintf(fid, '#\n');
% kappa = 0.01 = d phi/ dim_Z
% d_phi = w / (dim_X/2)
fprintf(fid, 'BoundaryCondition  3 loadTimeFunction 3 prescribedvalue %3.4d\n', 0.01*dim_X/2 * dim_Z);
fprintf(fid, '#\n');
fprintf(fid, 'BoundaryCondition  4 loadTimeFunction 4 prescribedvalue %3.4d\n', 0.01*dim_Y/2 * dim_Z);
fprintf(fid, '#\n');
fprintf(fid, '# load for indirect control\n');
fprintf(fid, 'NodalLoad 5 loadTimeFunction 5 Components 3 0. 0. -1.');
fprintf(fid, '\n');

% TIME FUNCTIONS

fprintf(fid, '#\n#\n#\n# TIME FUNCTIONS\n#\n#\n');

fprintf(fid, 'ConstantFunction 1 f(t) 1.0\n');
fprintf(fid, 'PiecewiseLinFunction 2 nPoints 3 t 3 -1.e6 0. 300. f(t) 3 0. 0. %3.3f\n', LOAD_eps);
fprintf(fid, 'PiecewiseLinFunction 3 nPoints 3 t 3 -1.e6 0. 300. f(t) 3 0. 0. %3.3f\n', LOAD_kappa1);
fprintf(fid, 'PiecewiseLinFunction 4 nPoints 3 t 3 -1.e6 0. 300. f(t) 3 0. 0. %3.3f\n', LOAD_kappa2);
fprintf(fid, 'ConstantFunction 5 f(t) 1.0\n');

% SETS

fprintf(fid, '#\n#\n#\n# SETS\n#\n#\n');

fprintf(fid, '#\n# concrete\n#\n');
fprintf(fid, 'set  1 elementranges {(%3d %3d)}\n', [1 n_elem_bricks]);

fprintf(fid, '#\n# large spiral\n#\n');
first = n_elem_bricks+1;
last = first+(length(SPIRAL(1).nr) - 2);
fprintf(fid, 'set  2 elementranges {(%3d %3d)}\n', [first last]);

fprintf(fid, '#\n# small spiral\n#\n');
first = last+1;
last = first;
for spiral = 2:nSpiral
    last = last + length(SPIRAL(spiral).nr) - 1;
end
last = last - 1;

fprintf(fid, 'set  3 elementranges {(%3d %3d)}\n', [first last]);

fprintf(fid, '#\n# vertical rebars\n#\n');
first = last+1;
last = first + nRebar * n_elem_Z - 1;
fprintf(fid, 'set  4 elementranges {(%3d %3d)}\n', [first last]);



% EXTRACTOR RECORD
fprintf(fid, '#\n# EXTRACTOR RECORD\n#\n');

fprintf(fid, '#%%BEGIN_CHECK%%\n');
fprintf(fid, '#TIME\n');
fprintf(fid, '#LOADLEVEL\n');
fprintf(fid, '#NODE number 1000000 dof 3 unknown d\n');
fprintf(fid, '#NODE number 1000001 dof 3 unknown d\n');
fprintf(fid, '#NODE number 1000002 dof 3 unknown d\n');
fprintf(fid, '#NODE number 1000003 dof 3 unknown d\n');

fprintf(fid, '#REACTION number 1000001 dof 3\n');
fprintf(fid, '#REACTION number 1000002 dof 3\n');
fprintf(fid, '#REACTION number 1000003 dof 3\n');

fprintf(fid, '#%%END_CHECK%%\n');

fclose(fid);


function  [choiceS, choiceSs, choiceL, choiceLs, choiceV, choiceVs] = rebar_dialog

screen_size = get(0,'screensize');

% d = dialog('Position',[600 500 250 250],'Name','Rebar selection');
d = dialog('Position',[screen_size(3)/2-200  screen_size(4)/2-125 450 500],'Name','Reinforcement specification');

% small

txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[10 400 275 40],...
    'String','Small spiral rebar diameter');

choiceS = '#3 (D10)';

popupS = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[50 400 100 25],...
    'String',{'#3 (D10)';'#4 (D13)';'#5 (D16)';'#6 (D19)';'#7 (D22)';'#8 (D25)';'#9 (D29)';'#10 (D32)'},...
    'Value',1,...
    'Callback',@popup_callbackS);


txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[10 350 210 40],...
    'String','Small spiral steel');

choiceSs = 'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%';

popupSs = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[50 350 325 25],...
    'String',{'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%';...
    'SD420 fy 497 MPa, fu = 723 MPa, eps = 10%';...
    'B500A fy 500 MPa, fu = 525 MPa, eps = 2.5%';...
    'B500B fy 500 MPa, fu = 540 MPa, eps = 5.0%';...
    'B500C fy 500 MPa, fu = 575 MPa, eps = 7.5%';...
    'fy 500 MPa, eps = infty'},...
    'Value',1,...
    'Callback',@popup_callbackSs);


% large

txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[10 280 275 40],...
    'String','Large spiral rebar diameter');

choiceL = '#4 (D13)';

popupL = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[50 280 100 25],...
    'String',{'#3 (D10)';'#4 (D13)';'#5 (D16)';'#6 (D19)';'#7 (D22)';'#8 (D25)';'#9 (D29)';'#10 (D32)'},...
    'Value',2,...
    'Callback',@popup_callbackL);

txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[10 230 210 40],...
    'String','Large spiral steel');

choiceLs = 'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%';

popupLs = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[50 230 325 25],...
    'String',{'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%';...
    'SD420 fy 497 MPa, fu = 723 MPa, eps = 10%';...
    'B500A fy 500 MPa, fu = 525 MPa, eps = 2.5%';...
    'B500B fy 500 MPa, fu = 540 MPa, eps = 5.0%';...
    'B500C fy 500 MPa, fu = 575 MPa, eps = 7.5%';...
    'fy 500 MPa, eps = infty'},...
    'Value',1,...
    'Callback',@popup_callbackLs);

% vertical

txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[10 160 250 40],...
    'String','Vertical rebar diameter');

choiceV = '#7 (D22)';

popupV = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[50 160 100 25],...
    'String',{'#3 (D10)';'#4 (D13)';'#5 (D16)';'#6 (D19)';'#7 (D22)';'#8 (D25)';'#9 (D29)';'#10 (D32)'},...
    'Value',5,...
    'Callback',@popup_callbackV);

txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[10 110 220 40],...
    'String','Vertical rebar steel');

choiceVs = 'SD420 fy 497 MPa, fu = 723 MPa, eps = 10%';

popupVs = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[50 110 325 25],...
    'String',{'SD420 fy 473 MPa, fu = 668 MPa, eps = 10%';...
    'SD420 fy 497 MPa, fu = 723 MPa, eps = 10%';...
    'B500A fy 500 MPa, fu = 525 MPa, eps = 2.5%';...
    'B500B fy 500 MPa, fu = 540 MPa, eps = 5.0%';...
    'B500C fy 500 MPa, fu = 575 MPa, eps = 7.5%';...
    'fy 500 MPa, eps = infty'},...
    'Value',2,...
    'Callback',@popup_callbackVs);


btn = uicontrol('Parent',d,...
    'Position',[89 20 70 25],...
    'String','Close',...
    'Callback','delete(gcf)');

% Wait for d to close before running to completion
uiwait(d);

    function popup_callbackS(popupS,event)
        idx = popupS.Value;
        popup_itemsS = popupS.String;
        
        choiceS = char(popup_itemsS(idx,:));
        
    end

    function popup_callbackSs(popupSs,event)
        idx = popupSs.Value;
        popup_itemsS = popupSs.String;
        
        choiceSs = char(popup_itemsSs(idx,:));
        
    end

    function popup_callbackL(popupL,event)
        idx = popupL.Value;
        popup_itemsL = popupL.String;
        
        choiceL = char(popup_itemsL(idx,:));
    end


    function popup_callbackLs(popupLs,event)
        idx = popupLs.Value;
        popup_itemsL = popupLs.String;
        
        choiceLs = char(popup_itemsLs(idx,:));
    end


    function popup_callbackV(popupV,event)
        idx = popupV.Value;
        popup_itemsV = popupV.String;
        
        choiceV = char(popup_itemsV(idx,:));
    end

    function popup_callbackVs(popupLs,event)
        idx = popupVs.Value;
        popup_itemsV = popupVs.String;
        
        choiceVs = char(popup_itemsVs(idx,:));
    end

end


function  [choiceBC, choiceLOAD] = loading_dialog

screen_size = get(0,'screensize');

d = dialog('Position',[screen_size(3)/2-125  screen_size(4)/2-125 300 180],'Name','Boundary conditions');

% BOUNDARY CONDITIONS

txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[20 120 210 40],...
    'String','Supports');

choiceBC = 'Periodic BC';

popupBC = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[25 120 250 25],...
    'String',{'Fully fixed';'Vertically fixed';'Periodic BC'},...
    'Value',3,...
    'Callback',@popup_callbackBC);


% LOADING

txt = uicontrol('Parent',d,...
    'Style','text',...
    'Position',[20 50 210 40],...
    'String','Loading');

choiceLOAD = 'Axial compression (direct)';

popupLOAD = uicontrol('Parent',d,...
    'Style','popup',...
    'Position',[25 50 250 25],...
    'String',{'Axial compression (direct)';'Eccentric compression (indirect)';'Pure bending (direct)';'Compression and bending (direct)'},...
    'Value',1,...
    'Callback',@popup_callbackLOAD);

btn = uicontrol('Parent',d,...
    'Position',[89 20 70 25],...
    'String','Close',...
    'Callback','delete(gcf)');


% Wait for d to close before running to completion
uiwait(d);

    function popup_callbackBC(popupBC,event)
        idx = popupBC.Value;
        popup_itemsBC = popupBC.String;
        
        choiceBC = char(popup_itemsBC(idx,:));
        
    end

    function popup_callbackLOAD(popupLOAD,event)
        idx = popupLOAD.Value;
        popup_itemsLOAD = popupLOAD.String;
        
        choiceLOAD = char(popup_itemsLOAD(idx,:));
    end

end






