# Example: 3-dimensional beam with T-section
#          Alternative definition of beam coordinate system
#
# --------------------------------------------------------------

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

# Data (N, mm)

  a =      1000;  % Length of beam
  b =        50;  % Width of cross section
  h =        50;  % Height of cross section
  t =         5;  % Thickness of cross section

  E   =   2.1E5;  % Young's modulus 
  ny  =     0.3;  % Poisson ratio
  rho = 7.85e-9;  % Mass density

  nel    =   20;  % Number of elements
  nmodes =   10;  % Number of modes

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

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

# Elements along beam

  geom   = mfs_beamsection("T", b, h, t);
  geom.Q = [0, 0, -1];

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

  idA = 1;
  idB = 1 + nel;

  nodes = struct("id",   {idA, idB},
                 "coor", {[0, 0, 0], [a, 0, 0]});

  [nodes, elem] = mfs_line(nodes, idA, idB, 
                           (idA + 1) : (idB - 1),
                           1 : nel, "b2", geom, mat);

# Graphic elements to visualize torsion

  ixC = 1 + floor(nel/2);
  idC = nodes(ixC).id;

  idB1 = idB + 1; idC1 = idB1 + 1;

  d = [0, 0, 0.1 * a];
  nodes(idB1).id = idB1; nodes(idB1).coor = nodes(2).coor + d;
  nodes(idC1).id = idC1; nodes(idC1).coor = nodes(ixC).coor + d;

  elem(nel + 1 : nel + 2) = ...
       struct("id",    {nel + 1, nel + 2},
              "type",  "g2",
              "nodes", {[idB, idB1], [idC, idC1]},
              "geom", [], "mat", []);

  model.nodes    = nodes;
  model.elements = elem;

# Rigid beams to connect graphic elements

  model.constraints.rigbdy = struct("noda", {idB, idC},
                                    "nodd", {idB1, idC1},
                                    "dofs", 1 : 3);

# Constraints

  model.constraints.prescribed = struct("id", idA, 
                                        "dofs", 1 : 6);

# Analysis
# --------

# Create and plot component

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

# Analysis

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

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

  fclose(fid);
