# Example: Thin-walled rectangular tube under torsion
#
# --------------------------------------------------------------------

# Data

  theta =   0.01;  % Torsion angle (in rad)
  a     =    100;  % Length in mm
  t     =      1;  % Thickness in mm
  E     = 210000;  % Young's modulus in MPa
  ny    =    0.3;  % Poisson's ratio

  b     = 0.5 * a;

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

# Output file

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

# Model definition

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

  % Nodal points at left end

  nodes( 1 :  4) = struct("id", {1, 2, 3, 4},
                          "coor", {[ 0, -a, -b], ...
                                   [ 0,  a, -b], ...
                                   [ 0,  a,  b], ...
                                   [ 0, -a,  b]});
  % Nodal points at right end

  nodes( 5 :  8) = struct("id", {101, 102, 103, 104},
                          "coor", {[ 10 * a, -a, -b], ...
                                   [ 10 * a,  a, -b], ...
                                   [ 10 * a,  a,  b], ...
                                   [ 10 * a, -a,  b]});

  % Nodal points for rigid bodies at left end

  nodes( 9 : 10) = struct("id", {1001, 1002},
                          "coor", {[ 0,  0, -b], ...
                                   [ 0,  0,  b]});

  % Nodal points for rigid bodies at right end

  nodes(11 : 13) = struct("id", {2000, 2001, 2002},
                          "coor", {[ 10 * a,  0,  0], ...
                                   [ 10 * a, -a,  0], ...
                                   [ 10 * a,  a,  0]});

  % Remaining nodal points

  nodes = mfs_linenodes(nodes, 1, 101, 11 : 10 : 91);
  nodes = mfs_linenodes(nodes, 2, 102, 12 : 10 : 92);
  nodes = mfs_linenodes(nodes, 3, 103, 13 : 10 : 93);
  nodes = mfs_linenodes(nodes, 4, 104, 14 : 10 : 94);

  tube.nodes = nodes;

  % Membrane elements

  ix = 1 : 10;
  for n = 1 : 4
      eltids = [n : 10 : 90 + n]';
      id4    = 5 - n;
      id1    = mod(id4, 4) + 1;
      nodids = [id1 : 10 : 90 + id1;
                id1 + 10 : 10 : 100 + id1;
                id4 + 10 : 10 : 100 + id4;
                id4 : 10 : 90 + id4]'; 
      quad(ix) = struct("id", num2cell(eltids),
                        "nodes", mat2cell(nodids, ones(10, 1)),
                        "type", "q4", "geom", geom,
                        "mat", mat);
      ix += 10; 
  endfor

  % Graphic elements

  graphic(1 :  4) = struct("id", {1001, 1002, 1003, 1004},
                           "nodes", {[1, 1001], [1001, 2], ...
                                     [3, 1002], [1002, 4]},
                           "type", "g2", "geom", [], "mat", []);  
  graphic(5 :  8) = struct("id", {2001, 2002, 2003, 2004},
                           "nodes", {[101, 2001], [2001, 104], ...
                                     [102, 2002], [2002, 103]},
                           "type", "g2", "geom", [], "mat", []);
  graphic(9 : 12) = struct("id", {3001, 3002, 3003, 3004},
                           "nodes", {[2000, 101], [2000, 102], ...
                                     [2000, 103], [2000, 104]},
                           "type", "g2", "geom", [], "mat", []);

  tube.elements = [quad, graphic];

  % Rigid bodies

  rigbdy(1) = struct("noda", 1001, "nodd", [  1,   2], "dofs", 1);
  rigbdy(2) = struct("noda", 1002, "nodd", [  3,   4], "dofs", 1);
  rigbdy(3) = struct("noda", 2001, "nodd", [101, 104], "dofs", 1);
  rigbdy(4) = struct("noda", 2002, "nodd", [102, 103], "dofs", 1);
  rigbdy(5) = struct("noda", 2000, "nodd",  101 : 104, 
                     "dofs", [2, 3]);

  tube.constraints.rigbdy = rigbdy;

  % Supports

  supports(1 : 4) = struct("id", {1, 2, 3, 4}, "dofs", [2, 3]);
  supports(5 : 6) = struct("id", {1001, 1002}, "dofs", 1 : 5);
  supports(7)     = struct("id", 2000, "dofs", 1 : 6);
  supports(8 : 9) = struct("id", {2001, 2002}, 
                           "dofs", [1 : 4, 6]);

  tube.constraints.prescribed = supports;

  % Prescribed rotation

  tube.loads.disp = struct("id", 2000, 
                           "data", theta * [0, 0, 0, 1]);

# Analysis

  cmp = mfs_new(fid, tube);
  mfs_export("tube.msh", "msh", cmp, "mesh");

  cmp = mfs_stiff(cmp);
  cmp = mfs_statresp(cmp);
  cmp = mfs_results(cmp, "statresp", "element");
  mfs_print(fid, cmp, "statresp", 
            "reac", "icsl", "disp", "stress");
  mfs_export("tube.pos", "msh", cmp, "statresp", "disp");

# Analytical solution

  geoma = mfs_beamsection("box", "thin", 2 * a, a, t);
  G   = E / (2 * (1 + ny));
  Mx  = G * geoma.IT * theta / (10 * a);
  tsx = Mx / (4 * a^2 * t);
  u   = Mx / (16 * a * t * G);

  fprintf(fid, "\nAnalytical Solution:\n");
  fprintf(fid, "   Mx = %10.3e, tauxy = %10.3e, ", Mx, tsx);
  fprintf(fid, "u = %10.3e\n", u);

  fclose(fid);
