# 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
  uz  =      10;  % Displacement amplitude

# Simulation paramters

  nofmod =   10;  % Number of normal modes
  fmax   =  400;  % Cut-off frequency

# Translation data

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

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

  geom = struct("t", t);

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

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

  BEAM.Motion = struct("type", "loads",
                       "name", "disp",
                       "data", [0, 0, uz]);

  BEAM.damping = struct("type", "ratios",
                        "data", {{{5, 0.05}, 0.02}});

# 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);

# Sets for response

  rpoints

  beam = mfs_newset(beam, "nset", "near", pntA, "Point_A");
  beam = mfs_newset(beam, "nset", "near", pntB, "Point_B");
  beam = mfs_newset(beam, "eset", "near", pntC, "Elem_C");
  beam = mfs_newset(beam, "eset", "near", pntD, "Elem_D");

  idA = mfs_getset(beam, "nset", "Point_A");
  idB = mfs_getset(beam, "nset", "Point_B");
  idC = mfs_getset(beam, "eset", "Elem_C");
  idD = mfs_getset(beam, "eset", "Elem_D");

  fprintf(fid, "  Point_A = %4.0d, Point_B = %4.0d\n",
          idA, idB);
  fprintf(fid, "  Elem_C  = %4.0d, Elem_D  = %4.0d\n",
          idC, idD);

# Save results

  save -binary modes.bin beam

  fclose(fid);
