# Example: Stiffened cylindrical box
#
# ----------------------------------------------------------------

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

# Data (N, mm)

  % Parameters

  nofmod  = 20;  % Number of normal modes
  opts.rbc = 8;  % Criterion to detect rigid body modes

  % Shell thicknesses

  ends  = struct("t", 5);  % Shell thickness of ends
  skin  = struct("t", 2);  % Shell thickness of skin
  floor = struct("t", 5);  % Shell thickness of floor

  % Cross section of frames and stringers

  S1.h = 20;    % Height of I-section
  S1.b = 20;    % Width of I-section
  S1.t =  2;    % Thickness of I-section

  frame    = mfs_beamsection("I", S1.b, S1.h, S1.t, S1.t, "rot");
  stringer = mfs_beamsection("I", S1.b, S1.h, S1.t, S1.t);

  % Cross section of stiffeners of front and end

  S2.h = 50;
  S2.b = 50;
  S2.t =  5;

  caps = mfs_beamsection("I", S2.b, S2.h, S2.t, S2.t);

  % Cross section of floor bars

  S3.h = 100;
  S3.b =  50;
  S3.t =   5;

  floor_bars = mfs_beamsection("box", "thin", S3.b, S3.h, S3.t);

  % Material data

  mat   = struct("type", "iso", "E", 70000, "ny", 0.34,
                 "rho", 2.7e-9);

# Model definition

  data = struct("type", "solid", "subtype", "3d");

  % Beam elements

  frame.v = [1, 0, 0];
  frame.P = [-0.5 * S1.h, 0];
  data.Frames = struct("type", "elements", "name", "b2",
                       "geom", frame, "mat", mat);

  stringer.Q = [5000, 0, 0];
  stringer.P = [0, -0.5 * S1.h];
  data.Stringers = struct("type", "elements", "name", "b2",
                          "geom", stringer, "mat", mat);

  caps.v = [-1, 0, 0];
  caps.P = [0, -0.5 * S2.h];
  data.Left_Stiffeners = struct("type", "elements", "name", "b2",
                                 "geom", caps, "mat", mat);
  caps.v = -caps.v;
  data.Right_Stiffeners = struct("type", "elements", "name", "b2",
                                 "geom", caps, "mat", mat);

  floor_bars.v = [0, 0, 1];
  floor_bars.P = [0, 0.5 * S3.h];
  data.Floor_Bars = struct("type", "elements", "name", "b2",
                           "geom", floor_bars, "mat", mat);

  % Shell elements

  data.Caps = struct("type", "elements", "name", "s4",
                       "geom", ends, "mat", mat);
  data.Skin = struct("type", "elements", "name", "s4",
                       "geom", skin, "mat", mat);
  data.Floor= struct("type", "elements", "name", "s4",
                       "geom", floor, "mat", mat);

# Analysis

  model = mfs_import(fid, "cylinder.msh", "msh", data);
  cyl   = mfs_new(fid, model);
  mfs_export("cylinder.axes", "msh", cyl, "mesh", "axes");

  cyl = mfs_stiff(cyl);
  cyl = mfs_mass(cyl);
  mfs_massproperties(fid, cyl);

  cyl = mfs_freevib(cyl, nofmod, opts);
  mfs_print(fid, cyl, "modes", "freq");
  mfs_export("modes.dsp", "msh", cyl, "modes", "disp");

  fclose(fid);
