# Example: Box Beam
#          Normal modes and error estimation
#
# -------------------------------------------------------------------

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

# Data (N, mm):

  E   =  210000;  % Young's modulus
  ny  =     0.3;  % Poisson's ratio
  rho = 7.85E-9;  % Mass density
  t   =       3;  % Thickness
  F   =       1;  % Load
  az  =    9810;  % Prescribed vertical acceleration (1 g)

  aK  =    1e-5;  % Parameter for Rayleigh-damping
  aM  =      20;  % Parameter for Rayleigh-damping

  nofmod =   20;  % Number of normal modes
  fmax   =  800;  % Maximum excitation frequency

# Translation data

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

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

  BEAM.Shell = struct("type", "elements", "name", "s4",
                      "geom", geom, "mat", mat);

  BEAM.Constraints = struct("type", "constraints",
                            "name", "prescribed",
                            "dofs", 1 : 6);

  BEAM.Load = struct("type", "loads",
                     "name", "point",
                     "data", [0, 0, F],
                     "lc", 1);

  BEAM.Motion = struct("type", "loads",
                       "name", "acce",
                       "data", [0, 0, az],
                       "lc", 2);

  BEAM.damping = struct("type", "Rayleigh",
                        "data", [aK, aM]);

  BEAM.A.type = "nodeset";
  BEAM.B.type = "nodeset";

# Import model

  model = mfs_import(fid, "beam.msh", "msh", BEAM);

# Create component and export axes

  beam = mfs_new(fid, model);
  mfs_export("beam.axes", "msh", beam, "mesh", "axes");

# Compute normal modes

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

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

# Estimate reduction error

  mfs_reductionerror(fid, beam, fmax);

# Compute modal effective masses

  mfs_meffmass(fid, beam);

# Elements for stress output

  beam = mfs_newset(beam, "eset", "near", [5,  40, 30], "C");
  beam = mfs_newset(beam, "eset", "near", [5, -40, 30], "D");

# Save results

  save -binary modes.bin fmax beam

  fclose(fid);
