# Example: Stiffened plate with attached point mass
#
# ----------------------------------------------------------------

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

# Data (N, mm)

  nofmod = 5;             % Number of normal modes

  geomm = struct("m", 1e-3);  % Point mass
  geoms = struct("t", 2);     % Shell thickness
  I     = struct("b", 10, "h", 10, "t", 1, "s", 1);

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

# Model definition

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

  geomb   = mfs_beamsection("I", I.b, I.h, I.t, I.s);
  geomb.v = [0, 0, 1];
  geomb.P = [0, 0.5 * I.h];

  data.Stiffeners = struct("type", "elements", "name", "b2",
                           "geom", geomb, "mat", mat);

  data.Plate = struct("type", "elements", "name", "s4",
                      "geom", geoms, "mat", mat);

  data.Point_Mass = struct("type", "elements", "name", "m1",
                           "geom", geomm);

  data.Supports = struct("type", "constraints",
                         "name", "prescribed", "dofs", 1 : 3);

  data.Rigfit = struct("type", "elements", "name", "g2");

  data.Rigfit_d = struct("type", "nodeset");  % dependent
  data.Rigfit_a = struct("type", "nodeset");  % autonomous

  data.Rigbdy = struct("type", "elements", "name", "g2");

  data.Rbar_1d = struct("type", "nodeset");
  data.Rbar_1a = struct("type", "nodeset");
  data.Rbar_2d = struct("type", "nodeset");
  data.Rbar_2a = struct("type", "nodeset");
  data.Rbar_3d = struct("type", "nodeset");
  data.Rbar_3a = struct("type", "nodeset");
  data.Rbar_4d = struct("type", "nodeset");
  data.Rbar_4a = struct("type", "nodeset");

  [model, nsets] = mfs_import(fid, "plate.msh", "msh", data);

# Linear Constraints

  model.constraints.rigfit = ...
          struct("nodd", nsets.Rigfit_d,
                 "dofa", {{nsets.Rigfit_a, 1 : 3}});

  Rbar_a = {nsets.Rbar_1a, nsets.Rbar_2a, ...
            nsets.Rbar_3a, nsets.Rbar_4a};
  Rbar_d = {nsets.Rbar_1d, nsets.Rbar_2d, ...
            nsets.Rbar_3d, nsets.Rbar_4d};
  model.constraints.rigbdy = ...
          struct("dofs", 1 : 6,
                 "noda", Rbar_a,
                 "nodd", Rbar_d);

# Analysis

  plate = mfs_new(fid, model);
  
  plate = mfs_stiff(plate);
  plate = mfs_mass(plate);
  mfs_massproperties(fid, plate);

  plate = mfs_freevib(plate, nofmod);

  mfs_print(fid, plate, "modes", "freq", "reac", "icsl");
  mfs_export("modes.dsp", "msh", plate, "modes", "disp");

  mem = mfs_meffmass(fid, plate);
  mpf = sqrt([mem.abs]);
  w   = 2 * pi * mfs_getresp(plate, "modes", "freq");
  R   = mpf .* w.^2;

  fprintf(fid, 
          "Reaction load resultants from effective masses\n\n");
  fdisp(fid, R);

  fclose(fid);
