# Example: 2-dimensional truss
#
# -------------------------------------------------------------------

# Data (N, mm)

  F      =       1;  % Unit load
  E      =  210000;  % Young's modulus in MPa
  rho    = 7.85E-9;  % Mass density in t/mm^3
  A      =     400;  % Cross section area in mm^2
  a      =    1000;  % Length in mm
  nofmod =      10;  % Number of normal modes
  fmin   =      10;  % Lowest excitation frequency
  fmax   =     400;  % Highest excitation frequency
  df     =       1;  % Frequency increment

# Output file

  fid = fopen("truss.res", "wt");

# Output device depending on operating system

  if (isunix)
     DEVICE = "svg";
  else
     DEVICE = "png";
  endif
  FORMAT = ["-d", DEVICE];
  EXT    = [".", DEVICE];

# Model definition
# ----------------

# Model type and subtype

  model = struct("type", "solid", "subtype", "2d");

# Nodal points

  nodes( 1).id =  1; nodes( 1).coor = [ 0, 0];
  nodes( 2).id =  2; nodes( 2).coor = [ a, 0];
  nodes( 3).id =  3; nodes( 3).coor = [ 0, a];
  nodes( 4).id =  4; nodes( 4).coor = [ a, a];
  nodes( 5).id =  5; nodes( 5).coor = [ 0, 2 * a];
  nodes( 6).id =  6; nodes( 6).coor = [ a, 2 * a];
  nodes( 7).id =  7; nodes( 7).coor = [ 2 * a, 2 * a];
  nodes( 8).id =  8; nodes( 8).coor = [ 3 * a, 2 * a];
  nodes( 9).id =  9; nodes( 9).coor = [ 0.5 * a, 2.5 * a];
  nodes(10).id = 10; nodes(10).coor = [ 1.5 * a, 2.5 * a];
  nodes(11).id = 11; nodes(11).coor = [ 2.5 * a, 2.5 * a];

  model.nodes = nodes;

# Geometry and material

  geom = struct("A", A);
  mat  = struct("type", "iso", "E", E, "rho", rho);

# Elements

  rod( 1).id = 1; rod( 1).nodes = [  1,  3];
  rod( 2).id = 2; rod( 2).nodes = [  1,  4];
  rod( 3).id = 3; rod( 3).nodes = [  2,  4];
  rod( 4).id = 4; rod( 4).nodes = [  3,  4];
  rod( 5).id = 5; rod( 5).nodes = [  3,  5];
  rod( 6).id = 6; rod( 6).nodes = [  4,  5];
  rod( 7).id = 7; rod( 7).nodes = [  4,  6];
  rod( 8).id = 8; rod( 8).nodes = [  5,  6];
  rod( 9).id = 9; rod( 9).nodes = [  6,  7];
  rod(10).id =10; rod(10).nodes = [  7,  8];
  rod(11).id =11; rod(11).nodes = [  5,  9];
  rod(12).id =12; rod(12).nodes = [  9,  6];
  rod(13).id =13; rod(13).nodes = [  6, 10];
  rod(14).id =14; rod(14).nodes = [ 10,  7];
  rod(15).id =15; rod(15).nodes = [  7, 11];
  rod(16).id =16; rod(16).nodes = [ 11,  8];
  rod(17).id =17; rod(17).nodes = [  9, 10];
  rod(18).id =18; rod(18).nodes = [ 10, 11];

  for n = 1 : length(rod)
      rod(n).type = "r2";
      rod(n).geom = geom; rod(n).mat  = mat;
  endfor

  model.elements = rod;

# Constraints

  model.constraints.prescribed = struct("id",   {1, 2},
                                        "dofs", [1, 2]);

# Loads

  model.loads.point = struct("id",   8, 
                             "data", {[F, 0], [0, F]},
                             "lc",   {1, 2});

# Damping

  model.damping = struct("type", "ratios",
                         "data", [0.05, 0.05, 0.02]);

# Analysis
# --------

# Create and plot component

  truss = mfs_new(fid, model);
  h = figure(1, "position", [100, 100, 750, 500],
                "paperposition", [0, 0, 15, 10]);
  subplot (1, 2, 1);
     mfs_plot(truss, "nodid", 1, "force", 1,
                     "loadcase", 1, "figure", h);
  subplot (1, 2, 2);
     mfs_plot(truss, "eltid", 1, "force", 1,
                     "loadcase", 2, "figure", h);
  print(["model", EXT], FORMAT, "-F:10");
  mfs_print(fid, truss, "loads", "point");

# Compute normal modes

  truss = mfs_stiff(truss);
  truss = mfs_mass(truss, "consistent");
  mfs_massproperties(fid, truss);

  truss = mfs_freevib(truss, nofmod);
  mfs_print(fid, truss, "modes", "freq");

# Estimate reduction error

  mfs_reductionerror(fid, truss, fmax);

# Plot some modes

  h = figure(2, "position", [200, 100, 500, 500],
                "paperposition", [0, 0, 15, 15]);
  subplot(2, 2, 1);
     mfs_plot(truss, "deform", 1, "modes", 1, "figure", h);
  subplot(2, 2, 2);
     mfs_plot(truss, "deform", 1, "modes", 2, "figure", h);
  subplot(2, 2, 3);
     mfs_plot(truss, "deform", 1, "modes", 3, "figure", h);
  subplot(2, 2, 4);
     mfs_plot(truss, "deform", 1, "modes", 4, "figure", h);
  print(["modes", EXT], FORMAT, "-F:10");

# Compute frequency response of both load cases

  f0 = fmin : df : fmax;

  [truss, f] = mfs_freqresp(truss, f0, "loadcase", 1);
  truss      = mfs_freqresp(truss, f, "nband", 0, "loadcase", 2);

# Get coefficients of first four normal modes

  Q1 = mfs_getresp(truss, "freqresp", "Q", 1 : 4, 1);
  Q2 = mfs_getresp(truss, "freqresp", "Q", 1 : 4, 2);

  figure(3, "position", [300, 100, 750, 750],
            "paperposition", [0, 0, 15, 15]);
  subplot(2, 1, 1);
     semilogy(f, abs(Q1));
     grid;
     legend('|Q_1 / F_x|','|Q_2 / F_x|','|Q_3 / F_x|',
            '|Q_4 / F_x|',
            'location', 'north');
  subplot(2, 1, 2);
     semilogy(f, abs(Q2));
     grid;
     legend('|Q_1 / F_y|','|Q_2 / F_y|','|Q_3 / F_y|',
            '|Q_4 / F_y|',
            'location', 'north'); 
     xlabel('f [Hz]');
  print(["Q", EXT], FORMAT, "-F:10");

# Get transfer functions of tip displacements

  rid = [8, 1; 8, 2];    % node id / dof

  U1 = mfs_getresp(truss, "freqresp", "disp", rid, 1); 
  U2 = mfs_getresp(truss, "freqresp", "disp", rid, 2); 

  figure(4, "position", [400, 100, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  semilogy(f, abs(U1(1, :)), "color", "green",
           f, abs(U1(2, :)), "color", "blue",
           f, abs(U2(1, :)), "color", "magenta",
           f, abs(U2(2, :)), "color", "red");
  grid;
  legend('|U_x / F_x|', '|U_y / F_x|', '|U_x / F_y|',
         '|U_y / F_y|',
         'location', 'north');
  xlabel('f [Hz]');
  ylabel('|U/F| [mm/N]');
  print(["U", EXT], FORMAT, "-F:10");

# Get transfer functions of normal forces in elements 1 to 3

  N1 = mfs_getresp(truss, "freqresp", "resultant", 1 : 3, 1);
  N2 = mfs_getresp(truss, "freqresp", "resultant", 1 : 3, 2);

  figure(5, "position", [500, 100, 750, 750],
            "paperposition", [0, 0, 15, 15]);
  subplot(2, 1, 1)
     semilogy(f, abs(N1{1}.N), "color", "green",
              f, abs(N1{2}.N), "color", "blue",
              f, abs(N1{3}.N), "color", "red");
     grid;
     legend('Rod 1', 'Rod 2', 'Rod 3',
            'location', 'southwest');
     ylabel('|N / F_x|');
  subplot(2, 1, 2)
     semilogy(f, abs(N2{1}.N), "color", "green",
              f, abs(N2{2}.N), "color", "blue",
              f, abs(N2{3}.N), "color", "red");
     grid;
     xlabel('f [Hz]');
     ylabel('|N / F_y|');
  print(["N", EXT], FORMAT, "-F:10");

# Save component for restart

  save -binary truss.bin truss

  fclose(fid);
