# Example: 2-dimensional frame
#
# ---------------------------------------------------------------

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

# Output device depending on operating system

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

# Data (N, mm)

  nmodes =   25; % Number of modes to be computed

  a      = 2500; % Length
  nela   =   20; % Number of elements per length

  E   =  210000; % Young's modulus
  rho = 7.85E-9; % Mass density

# Beam cross section data: Thin-walled square section

  b =        0.01 * a; % Length of edge of square
  t =         0.1 * b; % Thickness
  A =       4 * b * t; % Area
  I = 2 * b^3 * t / 3; % Moment of inertia

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

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

# Geometry and material

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

# Outer frame

  idA = 1;
  idB = idA + nela;
  idC = idB + nela;
  idD = idC + nela;
  idE = idD + nela;
  idF = idE + nela;
  idG = idF + nela;
  idH = idG + nela;
  idK = idH + nela;

  nodes(1).id = idA; nodes(1).coor = [0, 0];
  nodes(2).id = idD; nodes(2).coor = [0, 3 * a];
  nodes(3).id = idF; nodes(3).coor = [2 * a, 3 * a];
  nodes(4).id = idK; nodes(4).coor = [2 * a, 0];

  [nodes, elem1] = mfs_line(nodes, idA, idD, (idA + 1) : (idD - 1),
                            1 : (3 * nela), "b2", geom, mat);
  [nodes, elem2] = mfs_line(nodes, idD, idF, (idD + 1) : (idF - 1),
                            (3 * nela + 1) : (5 * nela),
                            "b2", geom, mat);
  [nodes, elem3] = mfs_line(nodes, idF, idK, (idF + 1) : (idK - 1),
                            (5 * nela + 1) : (8 * nela),
                            "b2", geom, mat);

# Cross beams BH and CG

  id1 = idK + 1;
  id2 = idK + 2 * nela - 1;
  [nodes, elem4] = mfs_line(nodes, idB, idH, id1 : id2,
                            (8 * nela + 1) : (10 * nela),
                            "b2", geom, mat);
  id1 = id2 + 1;
  id2 = id1 + 2 * nela - 2;
  [nodes, elem5] = mfs_line(nodes, idC, idG, id1 : id2,
                            (10 * nela + 1) : (12 * nela),
                            "b2", geom, mat);

  model.nodes    = nodes;
  model.elements = [elem1, elem2, elem3, elem4, elem5];

# Constraints

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

# Analysis
# --------

# Create and plot component

  frame = mfs_new(fid, model);
  mfs_plot(frame, "paperposition", [0, 0, 10, 10]);
  print(["frame", EXT], FORMAT);

# Analysis

  frame = mfs_stiff(frame);
  frame = mfs_mass(frame);

  frame = mfs_freevib(frame, nmodes);
  mfs_print(fid, frame, "modes", "freq", "reac");

# Plot some modes

  h = figure("position", [100, 100, 750, 1000], ...
             "paperposition", [0, 0, 16, 20]);
  subplot(3, 2, 1);
     mfs_plot(frame, "deform", 2, "modes", 1, "figure", h);
  subplot(3, 2, 2);
     mfs_plot(frame, "deform", 2, "modes", 2, "figure", h);
  subplot(3, 2, 3);
     mfs_plot(frame, "deform", 2, "modes", 3, "figure", h);
  subplot(3, 2, 4);
     mfs_plot(frame, "deform", 2, "modes", 4, "figure", h);
  subplot(3, 2, 5);
     mfs_plot(frame, "deform", 2, "modes", 6, "figure", h);
  subplot(3, 2, 6);
     mfs_plot(frame, "deform", 2, "modes", 16, "figure", h);
  print(["modes", EXT], FORMAT, "-F:12");

# Modal effective masses

  mfs_meffmass(fid, frame);

  fclose(fid);
